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LIST  OF  SYMBOLS 


Fortran 

Symbol 

Symbol 

Definition 

Units 

CAID 

c 

a 

starting  weight  of  ignition 
aid 

lb 

CPROP 

CHGWT 

c 

P 

starting  weight  of 
propellant 

lb 

g 

gravitational  constant 
(used  as  a conversion  factor) 

ft/sec^ 

H 

h 

c 

convective  heat  transfer 
coefficient 

Btu/in.^-°R 

SPACE2(3) 

h 

r 

radiative  heat  transfer 
coefficient 

Btu/in.^-°R 

TOTMOL 

“t 

total  moles  of  gas  in 
chamber 

lb -mol 

RMOL 

rate  of  change  of  total 
moles  of  gas  in  the  chamber 

mol /sec 

RATE 

r 

linear  burning  rate  of  the 
propellant 

in. /sec 

TICN 

'"is 

time  of  propellant  Ignition 

msec 

TMAX 

t 

pm 

time  of  maximum  pressure 

msec 

V(S) 

w 

a 

weight  of  ignition  aid  burned 

lb. 

Y(3) 

w 

al 

weight  of  ignition  aid 
combustion  products  in  tho 
chamber 

lb 

DY(5) 

• 

w 

a 

ignition  aid  mass  burning 

rate  (w  “AP”) 
a 

Ib/sec 

DY(3) 

• 

w , 
al 

rate  of  change  of  ignition 

aid  combustion  products  in 

the  system  (w  ."v  w /w  ) 

' al  a a n s 

Ib/sec 

Y(2) 

w 

i 

weight  of  igniter 
combustion  products  In  the 
chamber 

lb 
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Units 

DY(2) 

• 

w. 

1 

rate  of  change  of  igniter 

combustion  products 

(w.  = - w.  w /w  ) 

1 n s"^ 

Ib/sec 

DNN 

• 

w 

n 

mass  discharge  rate 
through  the  nozzle 

Ib/sec 

Y(6) 

w 

P 

weight  of  propellant  burned 

lb 

V(4) 

w , 
pi 

weight  of  propellant 
combustion  products  in  the 
chamber 

lb 

DY(6) 

* 

w 

p 

mass  burning  rate  of  the 
propellant 

Ib/see 

DV(4) 

• 

pi 

rate  of  change  of  propellant  Ib/sec 

combustion  products  in  the 

chamber  (vented  chamber  operation) 

(w  , o w - w , w /w  ) 
pi  p pl  n S'* 

Y(l) 

w 

V 

weight  of  air  in  tho  chamber 

lb 

W(l) 

« 

w 

r 

rate  of  change  of  air  in  the 

system  (w  « - w w /w  ) 

' ' r r n s 

Ib/soc 

w 

$ 

weight  of  gases  in  tho  chamber 

w =w  ♦w,  ♦w  , ♦ w , 
s r 1 al  pl 

lb 

tWSYS 

Vi 

s 

rate  of  change  of  the  weight 
of  gases  in  the  system 

Ib/sec 

OB 

X 

distance  burned  into  the 
grain 

in . 

Y(8-12) 

X 

n 

distance  burned  into  the 
grain  for  the  n-th  charge 
increment  (n  has  values  of  1-5 

in. 

) 

SPACUl(lO) 

\ 

effective  throat  area(sonic 
control  assumed) 

in.2 

SPACl-2(2) 

A 

w 

bomb  wall  surface  area 

in.  2 

CP 

C 

P 

heat  capacity  at  constant 
pressure 

8tu/lb. 
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Units 

CVP 

heat  capacity  at  constant 

Btu/lb  - °F 

V 

volume 

DCVP 

« 

rate  of  change  of  heat 

Btu/lb  - °F-sec 

V 

capacity  at  constant  volume 

CSTAR 

C* 

characteristic  discharge 
velocity 

ft/sec 

OD 

D 

initial  grain  diameter 

in. 

PD 

D 

P 

initial  perforation  diameter 

in. 

OOD 

D' 

instantaneous  grain  diameter 
(D'  = D - 2x) 

in. 

PPD.PRFD 

D* 

instantaneous  perforation 

in. 

P 

diameter 

(D'p  = Dp^2x) 

in.*" 

U 

U 

end  area  of  grain 

L* 

cv 

energy  of  gases  in  the 
chamber 

L* 

energy  from  combustion  of 

P 

the  propellant 

mi 

total  heat  loss  to  walls  of 

Btu 

L 

the  chamber 

ARP 

instantaneous  heat  loss  rate 

Btu/sec 

H 

\ 

average  heat  loss  rate 

Btu/ sec 

G 

convective  heat  loss  rate 

Btu/sec 

IffL 

Kt 

radiative  heat  loss  rate 

Btu/sec 

GL 

L 

initial  grain  length 

in. 

GGL 

L' 

instantaneous  grain  length 

in. 

• Variables  used  in  derivation  only 
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Units 

XMA 

«a 

molecular  weight  of  ignition 
aid  combustion  products 

NA 

XMI 

Mi 

molecular  weight  of  igniter 
combustion  products 

NA 

XMP 

Mp 

molecular  weight  of 

NA 

propellant  combustion 
products 

M 

Y* 

molecular  weight  of  air 

NA 

(defined  as  29) 

XMW 

Ms 

molecular  weight  of  gases 
in  the  system 

NA 

Y(7),P 

P 

pressure 

psia 

DP 

■ 

P 

rate  of  change  of  pressure 
with  respect  to  time 

psia/soc 

PMAX 

experimentally  measured 
maximum  pressux’O  in  the 
system 

psia 

P 

stagnation  pressure 

psia 

PT«EO 

•^OS 

theoretically  computed 
maximum  pressure  in  the 

psia 

system 

RSYS 

R 

s 

gas  constant  for  the  system 

in.-lb/ib°F 

DR 

• 

R 

s 

rate  of  change  of  gas 
constant  for  the  system 

in.-lb/lb-°: 

RU 

universal  gas  constant 

in.  -lb/ooi°l 

AAB.S 

instantaneous  propellant 

in.^ 

surface  area 

sec 
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Symbol 

Symbol 

Definition  Units 

AB 

■^tn 

instantaneous  surface  area 
of  nth  charge  increment 

. 2 
in. 

TSYS 

T 

s 

gas  temperature  in  the 
chamber  (computed  as  K;  note 
also  that  °R=°F  +459.69) 

°R 

TSYS 

■ T * 

St 

stagnation  temperature 

°F 

SPACE1(36) 

bomb  wall  surface  temperature 

^R 

TAIO,  T5 

"oa 

isochorie  adiabatic  flame 
temperature  of  the  ignition 
aid 

°R 

T6.TZD' 

isochorie  adiabatic  flame 
temperature  of  the  propellant 

BVOL 

empty  bomb  volume 

in.^ 

VSYS 

system  volume 

in.^ 

OVSVS 

• 

V 

s 

rate  of  change  of  sysiota 
volume^ 

if  a . w n-  nw* . ♦ +w  J / p 
s s sap 

in.'*/ sec 

WIO.W'OD 

w 

initial  propellant  web 

in. 

ALPHA 

a 

angle  used  in  csultiporforated  deg 
grain  surface  calculations 
(HPCISC)  defined  in  Appendix  C 

BETA 

B 

angle  used  in  ^*^SC,  defined 
in  Appendix  C 

deg 

(LAMl 

r 

ratio  of  specific  heats 

NA 

ETA 

n 

propellant  covolumc 

in.^/lb 

beta 

« 

n 

rate  of  change  of 
propellant  covolumc 

in.Vlb- 

TtiETA 

0 

angle  used  in  MPGSC,  defined 
in  Appendix  C 

deg 
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Symbol 

Symbol 

Definition 

Units 

RHO 

P 

solid  propellant  density 

lb/ in. 

CPGM2 

r 

a fimction  of  the  specific 
heat  ratio  defined  in 

NA 
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I.  INTRODUCTION 


Among  the  parameters  of  interest  to  the  interior  ballistician  are, 
the  burning  rate  and  chemical  energy  of  the  propellant  used  in 
propelling  charges.  These  parameters  must  be  knovm  to  predict  the 
performance  to  be  expected  from  a given  gun- ammunition  system. 

Reliable  data  on  the  chemical  energy  of  the  propellant  may  be  obtained 
from  thermodynamic  calculations  based  on  the  chemical  composition  and 
heats  of  formation  of  the  propellant  ingredients.^*^  Propellant  burning 
rates,  however,  cannot  be  predicted  with  similar  reliability  although 
various  efforts  at  predicting  burning  rates  from  chemical  compositions, 
5>6»7  thermochemistry®  and  chemical  kinetics  have  been  and  are  being 


ComeTt  M.  .4.J  FhDi  Theory  of  the  Interior  Ballistioa  of  GunSt  New 
Xorkj  John  Wiley  and  Sonst  Ina.y  1960. 

2 

Paul  G.  Baer  and  Jerome  M.  Frankie^  "The  Simulation  of  Interior 
Ballietia  Performanoe  of  Guns  by  Digital  Computer  Progranij " Ballistic 
Research  Laboratories  Report  No.  1186^  Decenber  1962^  AD  #299980. 

3 

A.  C.  Haukland  and  W.  M.  Burnett y "Sensitivity  of  Interior  Ballistic 
Performance  to  Propellant  Thermochemical  Parameters y " Proceedings  of 
the  Tri-Service  Gun  Propellant  Symposiwny  11-13  October  19? 2 y P%oat^nny 
Arsenaly  Dover y NJy  Vol  ly  Section  7.3. 

4 

Paul  G.  Baer,  Ingo  W.  Mayy  and  Jerome  M.  Frankley  "A  Comparison  of 
Several  Predictive  Approaches  in  Charge  Establishment  for  Large 
Caliber  Artillery  Systems y " Proceedings  of  the  11th  JANNAF  Combustion 
Meeting y Jet  Propulsion  Laboratory  Pasaaendy  CA.y  September  19?4y 
Vot  1;  C.P.I.A.  Publication  261y  December  19?4y  The  Chemical  Propulsion 
Information  Agency y Silver  Springy  Md. 

^William  H.  Tsohappaty  Lt.-Coly  Ord.  Dept.y  Text-Book  of  Ordnance  and 
Gunneryy  1st  Ed.y  New  Yorky  John  Wiley  and  SonSy  Inc,,  1917,  pp.  111-118. 

Q 

Albert  A.  Bennett,  PhD,  "Tables  of  Interior  Ballistics,"  Ordnance 
Department  Pamphlet  No.  2039,  April  1921. 

n 

D.  W.  Riefler  and  D.  J.  Lowery,  "Linear  Bum  Rates  of  Ball  Propellants 
Based  on  Closed  Bomb  Firings, " Ballistic ■ Research  Ldboratoriee  Contractor 
Renrt  No.  172,  August  1974.  (AD  #9217041) 

O 

Henn  Muraour,  "Cliemie  Physique  sur  la  Reaction  Entre  La  Temperature 
de  la  Explosion  d'  une  Poudre  et  sa  Vitesse  de  Combostion,"  Comp,  rend., 
Vol  187,  1928,  pp.  289-294. 
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pursued.  (See,  footnote,*).  Useful  burning  rate  data,  therefore, 

must  still  be  obtained  experimentally. 

Classically,  two  principal  methods  have  been  used  to  obtain  burning 
rate  data  for  propellants.  These  are  the  strand  burner  and  the  closed 
bomb.  Both  involve  the  burning  of  propellant  samples  in  steel  containers 
of  sufficient  strength  to  withstand  high  pressures.  In  a strand  burner, 
a single  sample  (strand)  of  propellant  is  burned  cigai'ette  fashion  under 
essentially  a constant  pressure.  In  practice,  this  is  achieved  by 
connecting  a ballast  volume  to  the  combustion  chamber  so  that  the  gases 
given  off  by  the  propellant  contribute  negligibly  to  overall  system 
pressure.  The  regression  of  the  burning  surface  is  measured  by  timing 
the  intervals  between  breaking  cf  fuse  wires  embedded  along  the  length 
of  the  propellant  sample.  Thus  each  experiment  yields  a single  average 
value  of  linear  burning  rate  for  a given  pressure.  To  obtain  burning 
rates  for  a propellant  over  a range  of  pressures,  a number  of  experi- 
ments are  required.  The  resultant  data  are  fitted  to  some  mathematical 
form  which  then  allows  computation  of  the  values  of  the  burning  rate 
at  intermediate  points.  The  method  is  straight-torward  (though  time 
consuming)  and  has  been  used  for  many  years,  espicially  by  the  rocket 
community. 

*Cuvv&ni  Army  efforts  on  . r r Hion  mode  I ing  are  oentered  in  the 
Fundo'ventalB  of  Comhmtx^,:  '"aek  of  the  Energotio  Materials  Teohnology 
Program  of  the  US  Army  Materiel  Development  and  Readiness  Command j 
Alexandria,  VA  22323. 

^Ref  1,  pp.  42-84  and  referenoe  therein, 

^^Sryae  L.  Craoford,  Jr.,  Clayton  liuggett,  and  J.  J,  MoBrady,  "The 
Mechanism  of  Burning  of  Double  Base  Propellants,"  J.  Phts,  and  Colloid 
Cham. , Vol  54,  1850,  pp.  854-862. 

Robert  G,  Parr  af\d  Bryoe  L.  Crawford,  Jr,,  "A  Physical  Theory  of 
But'ning  of  Double  Base  Rocket  Propellants,"  J.  Phys.  and  Colloid 
Chm. , Vol  54,  1850,  pp.  928-854. 

1 ■’ 

“0.  K.  Rtae  atvd  Robert  Ginell,  "The  Theoinj  of  Burmng  of  Double  Base 
Rocket  Powders,"  Phys.  and  Colloid  Chem.,  Vol  54,  1950,  pp.  885-81?. 

S.  Wilfong,  S.  S,  Pentwr,  and  P.  Daniols,  "An  Hi/pothesis  for 
Propellant  Burning,"  Phya.  and  Colloid  Chem.,  Vol  54,  1950, 
pp.  863-972. 

^'“R.  D.  Geckler,  "Mechanism  of  „ '■mbuation  of  Solid  Propellants," 

Selected  Combustion  Problems,  London,  Butterwoi'tha  Scientific 
Publications,  1954,  vp.  289-339. 

^^D.  B.  Spalding,  "The  Theory  of  Solid  and  Liquid  Propellants," 
embus tion  and  Flame,  Vol  4,  1960,  pp.  59-76. 


The  second  method  of  obtaining  propellant  burning  rate  data  involves 
the  closed  bomb.  In  the  closed  bomb,  a statistically  adequate  number 
of  propellant  grains  are  ignited  and  allowed  to  bum  in  a fixed  volmne 
under  the  pressure  exerted  by  the  propellant  combustion  gases.  The 
pressure  in  the  chamber  builds  up  rapidly  and  is  recorded  as  a function 
of  time.  From  the  pressure-time  data  it  is  possible  to  derive  linear 
burning  rat©  information  for  the  propellant  over  a range  of  pressures 
from  a single  experiment,  a marked  advantage  over  the  repeated  tests 
required  with  the  strand  burner,*  But  where  the  closed  bomb  technique 
gains  in  experimental  economy,  it  loses  in  the  complexity  of  the  data 
reduction  method.  The  problem  has  been  attacked  in  a variety  of  ways 
by  a number  of  authors. 17-20  The  earlier  papers  were  aimed  at 
providing  methods  for  computing  the  data  by  hand.  This  lead  to  the  use 
of  a variety  of  simplifying  assumptions  both  in  the  development  of  the 
theory  and  the  form  functions  used.  The  more  recent  papers 
were  aimed  at  computtsr  solutions  to  the  problem  and,  in  general,  provide 
a more  complete  treatment  of  the  phenomenon.  A brief  bibliography  of 
closed  bomb  bum  rate  reduction  methods  is  included  at  the  end  of  this 
report . 

^Alternately^  the  data  recorded  may  be  the  first  derivative  of  preaaure 
with  reapeot  to  time  va.  preaaure.  This  ia  generally  reduced  to  an 
average  value  of  dp/dt  (obtained  at  C, 250^0, 3750,0.500^  and  0.625  of 
the  maximum  preaaure)  which  when  compared  with  the  value  for  a 
reference  propellant  (obtained  under  idenPioal  oonditiona)  gives  an 
idea  of  the  burning  oharaateriaiia  to  be  expected  of  the  sample. 

The  "cpiiokneaa'*  and  '*i'elaPive  quiohxeaa"^^  valuea  ao  obtained  can  be 
quite  useful  for  oorrelationa  with  weapon  performance  aharaoteriatioa . 


7 /? 

Method  801,1.1,  ^'Quiahieaa  and  Force  Measurement  of  Propellant 
(Closed  Bomb  Method),*'  (Reviaed  SI  Oct  75),  Military  Standard  286B, 
Department  of  Defense,  Washington,  DC  SO 301. 

17 

C,  M.  Dickey,  "Determination  of  Burning  Charaateriatioa  of  Propellant, " 
E.f.  duPotit  de  Nemours  and  Company,  Bxploaivea  Department,  Bwnaide 
laboratory.  Memorandum  Report  No.  31  (File  BL-lSB-'lCl) , March  1943, 

1 s 

James  //,  Wiegand,  "A  MePixod  of  Calculation  of  the  Buining  Rate  of 
Powdeim  from  dP/dt  vs.  P Recoi'da  for  aloaed  Chambers,"  Ballistic 
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The  objective  of  the  present  work  was  to  generate  a comprehensive 
data  roductioii  method  providing  more  versatility  than  is  available  from 
any  of  the  previous  methods.  The  equations  which  were  developed  provide 
for  ;he  presence  of  igniter  and  ignition  aid,  as  well  as  the  ambient  air 
present  in  the  bomb.  A more  sophisticated  heat  loss  treatment  is 
included.  Other  unique  features  include  allowance  for  changes  in  the 
thermodynarai'  characteristics  of  the  combustion  products  and  of  web 
deviations  in  the  propellant  sample.  In  addition,  the  treatment  allows 
for  analysis  of  aata  from  vented  chamber  experiments.  The  theory  was 
implemented  in  a progrcim  written  for  an  available  laboratory  computer. 

An  attractive  feature  of  the  program  is  its  user  oriented  "question  and 
answer"  format  which  allows  the  user  to  readily  modify  input  variables 
and  to  make  decisions  on  the  various  data  reduction  options.  The  reduced 
data  are  available  for  examination  in  easily  interpretable  format  in  a 
matter  of  minutes. 

II.  OVERVIEW  AND  PROGRAM  CAPABILITIES 

The  rrogra;n  cupab’'itiv;s  are  outlined  in  the  following  section, 

In  the  form  listed  hore  it  requires  as  input  a suitable  data  file  of 
presrurt  vs  time  and  first  derivative  of  pressure  with  respect  to  time 
and  a "header"  or  set  of  parameters  describing  propellant  geometry  and 
physical  characteristics,  propellant  tiiermochemistry  and  experimental 
parameters.* 

A listing  of  the  input  requirements  is  contained  in  Appendix  A, 
program  operation  is  interactive  and  is  controlled  from  a terminal 
keyboard.  Two  modes  of  operation  are  possible,  standard  and  nonstandard. 
The  standard  option  provides  a "normal"  burning  rate  analysis  on  a more 
or  loss  routine  basis.  This  analysis  assumes  full  burning  of  all  of  the 
ignition  aid  before  ignition  of  tlie  propcMant,  simultaneous  ignition  of 
all  the  propellant  grains,  a fixed  set  of  dimen  ions  f^r  all  propellant 
grains,  constant  thermocljeroistry  for  the  combustion  products,  and 
linear  heat  loss  during  the  combustion  of  the  propellant  sample. 

The  second,  nonstandard,  mode  of  opsration  allows  for  operatijn  of 
the  program  with  a variety  of  options  for  special  jpplicaticna.  A 
sujiunary  comparison  of  the  "standard"  and  "nonstandard"  modes  of  airalysis 
is  presented  in  Table  i. 


^ Tho  data  file  io  obtaimd  by  opamtia^  en  the  i\xnge  data  with  a 
an.oothing  artd  diffeventiation  "SCHECK. " The  pyogrwit  and 

operation^  are  to  be  deaot'ibed  in  a fiitxa'e  x^port. 
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Burning  Surface  Capabilities  provided:  7 Allows  input  of  any  surface  area 

sphere;  cylinder;  and  function  in  tabular  form, 

single,  seven,  and  nineteen- 
perf orated  cylinders. 


Use  of  the  program  in  either  the  standard  or  nonstandard  mode  has 
been  simplified  by  allowing  the  use  of  propellant,  igniter,  and  experi- 
mental data  all  from  the  same  input  file.  Data  on  propellant  geometry, 
thermochemistry,  and  dimensions  are  all  recorded  as  "header"  information 
at  the  time  that  the  data  are  taken.  (This  is  achieved  using  a separate 
program,  CBDAP,  c^losed  bomb  data  Requisition  program.  The  program  will  be 
described  in  another  report. X The  result  is  that  conversion  of  the 
pressure-time  data  to  linear  burning  rates  may  be  performed  simply  without 
detailed  knowledge  of  program  operations.  In  performing  a standard 
analysis,  once  the  operator  has  entered  the  input  data  file,  he  is  freed 
from  providing  any  additional  input  except  for  specifying  whether  only 
the  central  portion  (from  10  to  80  percent  of  maximum  pressure)  or  all 
of  the  curve  is  to  be  plotted. 

If  changes  are  to  be  made  in  some  of  the  propellant  or  igniter 
data  or  if  special  handling  of  the  data  is  required,  the  nonstandard 
mode  of  analysis  is  employed.  In  this  mode,  the  pertinent  propellant 
and  igniter  data  in  the  data  file  are  displayed  and  the  opportunity  for 
changing  or  accepting  the  respective  value  is  presented.  In  addition, 
the  opportunity  of  selecting  each  of  the  options  in  Table  I is  presented 
to  the  user.  Once  these  decisions  have  been  made,  data  reduction  proceeds 
as  before. 

Program  output  consists  of:  (a)  a summary  sheet  or  header  describing 
the  sample,  experimental  parameters  and  data  on  the  maximum  pressure  and 
selected  values  of  the  derivative  of  pressure  with  respect  to  time 
(dP/dtj,  (b)  a tabular  listing  of  pressure,  time,  dP/dt,  burning  rate, 
web  and  surface  fraction  data,  (c)  superimposed  plots  of  pressure  (P) 
versus  timo  (t)  and  dP/dt  versus  t.  (d)  a plot  of  dP/dt  versus  reduced 
pressure  * ^^nd  (o)  a log-log  plot  of  bui-ning  rate  as  a function 

of  pvossux'o.  The  burning  rate  versus  pressure  plot  includes  a printout 
of  the  coefficient  and  the  exponent  in  the  equation 

r = AP*^  , (1) 

whoi'o:  r =2  linear  burning  rate, 

P = pressure, 

A « bux'ning  rate  coefficient,  and 
n = burning  t'ato  exponent. 


obtained  by  a least  squares  fit  of  the  data  over  the  desired  pressure 
range  as  well  as  statistical  data  on  the  "goodness  of  fit."  An  example 
of  program  output  is  given  in  Appendix  B. 
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III.  THEORY 


A.  SIMPLIFIED  DERIVATION 

The  basic  objective  of  the  analysis  is  to  derive  linear  burning 
rate  data  for  propellant  samples  from  pressure-time  histories  of  their 
buming  process  in  a closed  bomb.  Simply,  the  event  involves  the  con- 
version of  a solid  sample  composed  of  a large  number  of  grains  of  a 
given  geometry  and  size  to  a gas  having  a given  amount  of  energy.  Since 
the  vessel  is  closed,  the  products  may  not  escape,  the  pressure  builds 
up  and  the  propellant  sample  bums  in  an  environment  of  a continuously 
changing  pressure. 

A variety  of  factors  influence  the  conversion  rate  of  the  solid 
sample  to  gas.  Those  of  primary  interest  are  the  propellant  surface 
area,  the  pressure,  and  the  propellant  chemical  composition.  For  all 
propellants,  the  conversion  rate  of  solids  to  gas  (i.e.,  the  mass  burning 
rate)  is  directly  proportional  to  surface.  For  all  gun  propellants  (and 
many  others)  the  rate  of  regression  of  the  propellant  surface  (the  so- 
called  linear  buming  rate,  r)  is  directly  proportional  to  pressure. 

As  a general  rule,  the  linear  burning  rates  of  propellants  at  given 
pressures  are  a function  of  the  energy  content  of  the  composition. 

To  describe  the  process  it  is  necessary  to  describe  the  gas  pi'oduced, 
the  unburned  propellant  and  the  energy  balance  for  the  process  as  a 
function  of  time.  The  derivation  which  follows  is  limited  to  a single 
propellant  of  constant  thermochemistry.  Eliminating  complicating  factors 
such  as  the  presence  of  the  igniter,  ignition  aid,  etc.,  allows  the 
generation  of  a simple  instructive  sot  of  equations  demonstrating  the 
logic  used.  The  actual  equations  used  in  the  program  are  discussed  in 
a later  section  and  their  derivation  is  given  in  the  appendix. 

(1)  Equation  of  State  of  Gas 

The  combustion  products  may  bo  described  using  the  following  equation 
of  state; 


whore : 

» system  volume 

Wp  » weight  of  propellant  burned 
R^  “ gas  constant  for  the  system 
T^  o gas  temperature  in  the  chamber 


(2) 
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The  equation  is  formally  analogous  to  the  familiar  Ideal  Gas  Equation. 
The  difference  is  in  the  definition  of  the  system  volume  term  (V  ) which 
is  defined  as  the  chamber  volume  modified  to  reflect  the  presence  of 
unburned  propellant  and  the  covolume  correction,  Equation  (3) 


V W 


p p 


where : 


= chamber  volume 

Cp  = initial  weight  of  propellant 
p = solid  propellant  density 
n = propellant  covolume. 


(3) 


It  should  be  noted  that  the  mixture  of  gases  making  up  the  propell- 
ant combustion  products  is  treated  as  if  it  were  a single  gaseous 
species  having  specific  properties  of  molecular  weight,  heat  capacity 
and  covolume.  These  properties,  of  course,  are  determined  by  the  nature 
and  stoichiometry  of  the  combustion  products. 

Once  the  appropriate  substitutions  are  made,  the  Equation  of  State 
becomes : 


P 


( 1 - n) 

p 


° w R T 
p s s 


(2)  Energy  Balance  Equation 


(4) 


The  energy  fx*om  combustion  of  the  propellant  sample  is  partitioned 
between  the  internal  energy  of  the  product  gases  and  heat  loss  to  the 
chamber.  The  Energy  Balance  Equation  may  bo  written  as; 


(5) 


whore;  E^^  « energy  of  gases  in  the  chamber 

0 a energy  from  combustion  of  the  propellant,  and 
P 

a heat  loss  to  walls  of  the  chamber 
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The  equation  may  be  rewritten  as: 


C.,w  T 
V p s 


SVop  ■ “l 


where:  Cy  ^ heat  capacity  at  constant  volume  and 


Tqp  = isochoric  flame  temperature  of  propellant. 
R 

s 

Pq^  = Theoretically  computed  maximum  pressure 
= Experimentally  measured  maximum  pressure 


The  temperature  of  the  combustion  gases,  T^,  is  less  than  the  isochoric 


adiabatic  flame  temperature  computed  for  the  propellant  composition,  T 
to  heat  losses  to  the  walls  of  the  chamber.  The  heat  capacity  at 
constant  volume,  Cy,  is  an  average  property  between  T^  and  T^^  for  the 

mixture  of  gases  making  up  the  combustion  products  of  the  formulation. 
It  is  defined  per  unit  weight,  rather  than  per  mole. 


Op’ 


due 


(3)  Rate  of  Conversion  of  Solid  to  Gas. 


To  obtain  the  equation  for  the  rate  of  conversion  of  the  soha 
propellant  to  gaseous  combustion  products.  Equations  (4)  and  (6)  arc 
differentiated.  Differentiation  of  Equation  (4)  holding  V.  , c , p,  n and 
R constant,  yields:  ^ 


P 


dP 

dt 


^ P ( i - n) 

p 


dW  a 


R w 
s p 


dT 
s 

dt 


R T 

s s 


(7) 


The  rate  of  conversion  of  solid  to  gas  is  given  by  dwp/dt,  the  rate  of 
formation  of  propellant  combustion  products.  This  is  the  term  we  ai'c 
seeking  to  evaluate  in  terms  of  experimental  parameters.  To  do  this, 
it  is  necessary  to  define  the  rate  of  chojigo  of  system  temperature 
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(dT  /dt) . This  is  done  by  differenting  the  Energy  Equation.  Differenti- 
atiSn  of  Equation  (6)  holding  and  Tqp  constant  and  rearranging,  yields 


(T 


2R. 


- T ) 


W 


dw 

__E 

dt 


!!k 


at  this  point,  the  following  relationship  is  introduced: 


(8) 


R3  = Cy  Cy  - 1) 


where : 


Y = ratio  of  heat  capacities 

The  relationship  is  strictly  true  for  ideal  gases  but  is  conunonly 
used  in  describing  real  systems.  It  allows  recasting  Equation  (8)  in 
the  following  form: 


(T 


dt  RIvT 

S P 


(9) 


Substituting  the  right  hand  side  of  Equation  (9)  into  the  differentiated 
Equation  of  State,  Equation  (7);  yields: 


V.  - c ^ w (1  - n) 
b ^ P _ 


dP 

dt 


Vs  - f 0 - n) 


dw 

ik 


♦ R W 
s p 


(T 


2e. 


dt 


(10) 


The  equation  may  now  bo  solved  for  dw  /dt,  giving  the  rate  of 
foiination  of  propellant  combustion  products  in  terms  of  experimental 
data  ( P,  dP/dt)  and  a series  of  constants  (Vj^,  c^,  p,  n,  R^, 

and  y) • The  resulting  equation  is : 


I 
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(11) 


dw 

dt 


c , w (1  - n)  dp  ^ (v-1) 

P P dt 


- Ml  - n) 

^ P 


(4)  Linear  Burning  Rate. 

The  linear  burning  rate,  r,  is  the  rate  of  regression  of  the 
propellant  surface  (dx/dt).  To  compute  it  one  begins  with  considering 
the  volume  element  burned  through  during  an  infinitesimal  time  interv'al . 
The  following  equation  applies: 


dWp  = p dx  (12) 

dt  dt 


where: 


» the  surface  area  of  the  propellant  at  any  time  t. 

dx  = rate  of  regression  of  the  propellant  surface,  equals  r, 

3T  the  linear  burning  rate. 

Equation  (12)  defines  the  rate  of  formation  of  propellant  combustion 
products  in  geometric  terms.  If  the  propellant  is  composed  of  a number 
of  identical  grains,  the  propellant  surface  area  may  bo  computed  using 
a variety  of  "form  function"  equations.  Equations  have  been  developed 
for  spherical,  cylindrical  and  perforated  cylindrical  (1,  7,  and  19 
perforations)  grain  geometries.  The  generalized  equation  is; 


S,  . f Cx) 


(13) 


where: 


X « the  distance  burned  into  the  grain. 

When  X equals  zero,  the  surface  area  is  the  initial  surface  area  of  the 
propellant  grain.  As  x increases  positively,  the  surface  area  of  the 
gi'ain  changes  characteristically  for  each  grain  typo.  Form  functions 
for  all  of  the  grain  typos  mentioned  above  have  been  included  in 

Appendix  C. 
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Examining  Equation  (12),  it  is  evident  that  the  objective  of 
computing  the  linear  burning  rate  of  the  propellant  from  the  closed  bomb 
pressure-time  data  has  been  attained.  The  term  dWg/dt  is  defined  by 
Equation  (11),  and  the  surface  area,  S , is  defined  by  the  Form  Function 
Equation,  Equation  (13) . This  completes  the  simplified  derivation  intended 
for  inclusion  in  the  text. 

B.  Theory  Used  in  The  Program. 

Both  the  event  described  and  the  equations  describing  it  are 
considerably  more  complex  than  just  described.  In  every  experiment  an 
igniter  is  used  to  start  the  propellant  burning  and  in  many  an  ignit- 
ion aid  is  used  as  well.  Further,  the  volume  in  the  bomb  not  occupied 
by  the  propellant  at  the  start  of  the  experiment  is  occupied  by  ambient 
air.  It  is,  therefore,  evident  that  one  is  not  dealing  with  the  single 
component  system  described  earlier.  The  Equation  of  State  used  in  the 
program  includes  all  components.  The  Energy  Equation  is  also  more 
complex,  since  the  thermodynamics  of  the  combustion  gases  are  allowed 
to  change  and,  in  addition,  allowance  is  made  not  only  for  heat  loss 
from  the  system  but  mass  loss  as  well  (vented  chamber  operation). 

The  Mass  Burning  Rate  Equation  derived  for  the  system  reflects  these 
complexities.  See  Equation  (27),  Appendix  D.  Since  the  treatment  allows 
for  deviations  in  the  ignition  of  the  propellant  charge  as  well  as  web 
deviations,  computing  the  surface  area  at  any  instant  is  also  slightly 
more  complicated.  Finally,  the  instantaneous  Heat  Loss  Terra  (H^^)  is 

evaluvited  in  the  pi'ogram  in  one  of  two  ways,  either  as  an  avoi'agc  value, 
constant  throughout  the  but'n,or  as  a variable  defined  by  the  radiative 
and  convective  heat  loss  elements.  The  objective  of  this  se-tion  is  to 
provide  some  explanatory  comments  on  several  equations  used  in 
the  program. 

(1)  Rato  of  Conversion  of  Solid  to  Gas. 

For  the  purpose  of  discussion  the  Mass  Burning  Rate  Equation  derived 
in  Appendix  D has  been  I'egroupcd  and  the  numerator  divided  into  a series 
of  texms  A through  E.  This  is  the  form  in  which  it  appears  below. 


dw 

P 


S 


A + B^C4-D>E 


(14) 


where ; 


=»  universal  gas  constant 

oolocuiar  weight  of  propellant  combustion  products 


a molecular  weight  of  gas  in  the  system 
A = V dP 

* at 


B a (Y-i) 


where : 


C = 


yR,T,.Pn 


dw 

n 

■It 


w 1 a weight  of  propellant  combustion  products  in  the  chamber 
^ (as  opposed  to  the  total  weight  of  propellant  burned) 


dw^  = gas  discharge  rate  through  the  nozzle 

"at 


D a W^T 
s s 


« Cp  it 


- Pw„  dn 
> 

3t 


c • 


P(n-  1) 

p 


-Oi  Hi  * ».  T,  (1  J>1 

a s J 


d w 


al 


It 


where : 


= molecular  weight  of  ignition  aid  combustion  products 
“al  * of  ignition  aid  combustion  products  in  the  chamber. 


Several  of  the  terms  are  as'cociatcd  with  exercising  program  options 
previously  described  (Table  I).  ITie  functions  of  Terms  A through  E are 
listed  in  Table  II. 
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Table  II.  Functions  of  Terms  A through  E in  Equation  (14). 

A System  volume  term 

B Heat  Loss. 

C Mass  Loss. 

D Variable  thermochemistry. 

E Contribution  of  simultaneously  burning  ignition  aid. 

The  System  Volume  Term  (A)  is  necessary  for  all  computations.  The 
Heat  Loss  Term  (B)  is  included  as  required  in  the  analysis.  The  term 
may  be  evaluated  simply  (standard  option,  Table  1)  or  comprehensively 
(nonstandard  option).  This  will  be  discussed  more  fully  in  section  B(3) . 
In  the  case  of  vented  chamber  operation  or  in  computing  burn  rates  from 
artificially  generated  pressure  time  data,  term  B can  be  zero.  The 
Mass  Discharge  Term,  (C)  is  used  in  analyzing  vented  chamber  experiments. 
Term  (D)»Variable  Thermochemistry,  is  important  in  analyzing  low  pressure 
closed  bomb  data  since  the  thox-mochemistry  of  the  combustion  products 
changes  significantly  with  pressure  at  low  pressures.  The  inclusion  of 
the  Ignition  Aid  Burning  Tenn  (E)  becomes  important  when  descri>^ing 
situations  involving  simultaneous  bunxing  of  px'opollant  and  ignitor. 

This  is  the  case  more  often  than  not,  though  the  decision  on  the  over- 
lap of  ignition  aid  and  px'opollant  buxixing  is  made  on  the  basis  of 
experience  by  the  pi'ogram  opei-ator. 

In  essence,  Equation  (14)  is  Equation  (11)  appropriately  modified 
to  vofloct  the  complexities  of  the  o-xpci-iment.  This  nii;,y  bo  x'oadily 
demonstrated  by  imposing  the  same  assumptions  on  Equation  (14)  as  were 
used  in  deriving  the  simplified  Mass  Burning  Rate  Equation  (Equation  U). 
Assuming  a single  px-opeUant  (no  igniter  or  ignition  aid),  con.stant 
thei-mochemistry,  closed  bomb  operation  (no  mass  loss  through  a nozzle) 
and  the  absence  of  air  from  the  chaadjcr,  the  following  terms  in  Equation 
(14)  may  be  eliminated: 

Term  C.  Under  closed  boad)  conditions  dw  /dt=0 

n 

Tcim  0.  For  const&nt  thermochemistf)’  dn/dtsQ 

Term  E.  In  the  absotxce  of  an  ignition  dw  ,/dt°0 
aid 


R T,  (1_  - 1 )|.  For  a single  component 

^ Np  system  (i/M^  - 1/M^)  = 0 
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Elimination  of  the  terms  above  reduces  Equation  ^14)  to  Equation 

(11) .  The  methods  used  in  deriving  ti;  - two  equations  were  the  same. 

The  derivation  of  Equation  (14)  paralleling  the  approach  used  in  the  • 
text  for  Equation  (11)  is  given  in  Appendix  D. 

(2)  Linear  Burning  Rate 

The  linear  burning  rate  in  the  program  is  computed  by  Equation 

(12) .  Differences  arise,  however,  in  the  computation  of  the 
instantaneous  burning  surface  area  S . Two  of  the  program  options 
treat  ignition  deviations  in  the  propellant  charge  and  web  deviation  in 
the  propellant  grains.  Both  options  influence  the  burning  surface  area. 

If  an  ignition  deviation  takes  place,  parts  of  the  propellant  charge 
begin  to  bum  before  others.  In  this  case  the  simple  computation  of 
total  burning  surface  area  as  a function  of  distance  burned  is  not 
appropriate.  What  is  done  is  to  proportion  the  propellant  charge  into 
five  parts  (two  each  of  10  5 jrcent,  two  each  of  20  percent  and  one 
of  40  percent)  and  to  allow  the  ignition  of  the  charge  increments  to 
differ  by  some  arbitrary  :ime  input  by  the  operator.  The  distances 
burned  into  the  surface  of  each  portion  of  the  charge  are  carried 
separately.  Under  these  conditions  the  burning  rurface  is 
computed  according  .0  the  following  equation: 

■ ^ ^ ! S (15) 

^ 1 ’•n 


where : 

S = surface  area  of  the  nth  propellant  charge  increment 

% 


= distance  burned  into  the  surface  of  the  nth  propellant 
charge  increment. 


Of  course,  for  simultaneous  ignition  of  the  propellant  charge, 
Equation  (15)  reduces  tc  Equation  (13).  It  must  be  emphasized,  however, 
that  one  has  no  a priori  knowledge  of  the  ignition  deviation  time  of  the 
propellant  charge;  so  this  treatment  should  be  viowod  with  caution. 

The  web  deviation  is  handled  analogously;  In  this  case,  however, 
web  deviation  values  may  be  obtained  from  actUul  measurements. 
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(3)  Heat  Loss 

(a)  Standard  Option.  Evaluation  of  td'*  Heat  Loss  Term  (H. ) is  done  in 
either  of  tv;o  ways.  In  the  standard  option  it  is  some  suitable  average 
value,  constant  throughout  the  analysis.  The  following  equation  applies: 


H,  = H,  = C,,V, 
L L Vs 

T 


('■orV 

(‘pni''ig) 


where : 

H. 

L 


= average  heat  loss  rate 


(16) 


Pq  = theoretically  con5>uted  maximum  pressure 
^ (adiabatic  conditions,  contributions  from  propellant, 
igniter,  ignition  aid  and  air  in  system) . 

tp^  = time  of  maximum  pressure 

t.  a time  of  ignition 


The  total  heat  loss  is  the  difference  betwee:'  the  adiabatically 
computed  internal  energy  of  the  system  and  the  internal  energy  computed 
from  the  maximum  pressure  observe^.  The  total  heat  loss  is  converted 

into  the  average  Heat  Loss  Rato  (H, ) by  dividing  by  the  burning  time 

interval  (t^  -t,  ) . 

pm  ig^ 

(b)  Nonstandard  Option 

In  the  nonstandard  of'ion,  heat  loss  is  analyzed  into  its  convective 
and  radiative  components.  It  is  assumed  that  during  the  time  that  the 
propellant  burns  the  gas  generation  results  in  convective  heat  transfer 
to  the  chPinbor  walls.  ITiis  is,  of  course,  accompanied  by  radiative  heat 
losses  as  well.  After  the  propellant  is  consumed  it  is  assumed  that  the 
convective  heat  lo^s  becomes  insignificant  relative  to  the  radiative 
heat  loss.  These  assumptions  allow  the  following  treatment. 

(i)  Radiative  heat  loss  coefficient.  The  temperature  of  the  gases  in  the 
system  is  given  by; 


T 

s 


a p V 
m s 


The  radiative  heat  loss  rate  is  given  by: 
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dt 
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where; 

» radiative  heat  loss  rate, 

Once  a matched  array  of  T , H.  data  are  generated,  they  may  be  fitted  to 

S LX 

the  following  relationship: 
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H, 


Lr 

“T“ 
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w 


where : 


(18) 


h a radiative  heat  transfer  coefficient 
r 

e bomb  surface  area 
= bomb  wall  temperature. 


In  carrying  out  these  computations,  T may  bo  neglected  since  it  is 

indeed  insignificant  (T^  »T^  ) in  the  radiant  heat  transfer.  The  bomb 

surface  area  (A^  ) is  one  of  the  program  inputs  (alternately  a default 

value  of  18.1  in^  is  available  in  the  program)  and,  therefore,  the 
evaluation  of  the  Radiant  Heat  Transfer  Coefficient  (h  ) is  accomplished. 
Throughout  the  analysis,  h is  a constant.  ^ 


(ii)  Convective  Heat  Transfer  Coefficient.  The  Convective  Heat  Ti'ansfor 
Coefficient  (h^)  is  assumed  to  bo  a function  of  the  instantaneous  mass 
flow  as  in  the  case  of  heat  transfer  in  a pipe.  The  relationship  is 
given  as: 


(dwj 

h ’F 
c c dt 


0.8 


(19) 


where : 


h * convective  heat  transfer  coefficient 
F a proportionality  constant 
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The  value  of  h is  computed  at  every  point  in  the  analysis  using  the 
instantaneous  mass  generation  rate.  The  value  of  the  proportionality 
constant  F is  obtained  using  an  approximation  technique  in  which  an 
approximate  value  of  h is  calculated  from  the  average  mass  generation 
and  heat  loss  rates  ana  this  approximate  value  is  refined  using  the 
value  of  h^  and  increments  of  AP/At. 

(iii)  Heat  Loss  Rate,  nonstandard  option.  To  compute  the  Heat  Loss  Rate 
(H^)  at  any  given  point  in  the  analysis  the  following  equations  are  used: 


H.  = h^  A (T  - T ) 
Lc  c w s w"^ 


where:  H.  = convective  heat  loss  rate 
Lc 


H, 


Lr 


4 4 

h , A (T  - T,  ) 
r w ^ s w 


(20) 


(21) 


H 


L 


(22) 


The  initial  value  of  the  Wall  Temperature  (T  ) is  assumed  to  be  450°K 
(this  may  be  changed  by  the  operator)  and  T^”is  continuously  adjusted 
throughout  the  analysis  as  is  the  Systems  Temperature  (T^).  The 
relationships  governing  heat  conduction  to  the  wall  are  given  in 
Appendix  E. 


IV.  PROGRAM  STRUCTURE  AND  OPERATION 
A.  Standard  Analysis 

The  program  consists  of  a main  program  and  a number  of  subroutines 
which  handle  reading  in  of  data,  unit  conversions , heat  loss  computations, 
burning  rate  calculations,  and  printing  and  plotting  of  the  output  data. 
The  program  structure  is  outlined  in  Figure  1.  The  diagram  indicates  the 
relationship  of  the  subroutines  in  the  program.  Capsule  summai'ies  of  the 
functions  of  the  subroutines  are  contained  in  Appendix  F. 

Program  operation  is  most  easily  described  by  following  a standard 
analysis  sequence  step  by  step.  The  options  may  then  be  seen  as  pertur- 
bations on  the  normal  mode  of  analysis.  A flow  diagram  of  the  program 
appears  in  Figure  2.  The  chart  has  been  arranged  to  show  the  operations 
involved  in  a standard  analysis  in  sequential  form.  The  options  possible, 
heat  loss,  mass  loss,  variable  thermochemistry,  etc.,  are  shown  offset 
from  the  main  sequence  of  operations.  The  optional  sections  are 
marked  with  arterisks. 
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PTOgram  Structure  CCBRED) . Interrelation  of  Subroutines . 


1 


CALL  CBRED 


Figure  2-1.  Closed  Bomb  Bum  Rate  Program  (CBRED). 
Generalized  Plow  Scheme. 
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Figure  2-3.  Closed  Bomb  Burn  Rato  Program  (CBRED) . 
Generalized  Flow  Scheme. 
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II  [If  III  Ji  ii|i^ 


1.  Start  up  (Blocks  1,  2 and  3).  Once  the  program  (CBRED)  is  called, 
it  signs  in  and  requests  the  file  name  and  extension  of  the  data 
file  to  be  analyzed.  On  input  of  the  file  name  and  extension  the  data 
file  is  read  in.  The  input  data  consists  of  matched  arrays  of  pressure- 
time and  time  derivative  of  pressure-time  data  as  well  as  a data  set 
consisting  of  experimental,  thermodynamic  and  grain  geometry  data 
(see  Appendix  A) . The  data  is  stored  on  magnetic  tape  and  is  handled 
as  a single  unit.  (At  the  time  the  experiments  run,  the  operator  records 
in  the  da ''a  file  all  the  parameters  required  for  data  reduction.  This 
is  done  using  a separate  program,  CBDAP.  The  pressure-time  data 
obtained  is  differentiated  by  a second  program,  SCHECK.  The  "header" 
information  is  carried  through  to  the  new  data  file.  Opportunities 
for  updating  the  "header"  information  are  provided  by  both  programs) . 
These  operations  take  place  in  the  main  program  and  subroutine  ACQUIRE. 

2.  Mode  of  Analysis  (Block  4).  The  question  is  posed  to  the  operator 
whether  the  analysis  to  be  performed  will  be  in  the  standard  or 
nonstandard  mode.  Once  the  choice  for  the  standard  analysis  mode  has 
been  made,  the  program  goes  into  execution  without  further  input  from 
the  operator.  See  subroutine  ACQUIRE  also. 

3.  Initialization  (Block  5).  At  this  point  a variety  of  operations 
necessary  to  initialize  the  problem  are  performed.  The  input  parameters 
are  converted  to  a consistent  set  of  units.  At  present  a pound,  foot 
(inch),  second,  BTU  system  is  employed,  but  gas  temperatures  are  in 
degrees  Kelvin.  A number  of  constants  of  multiple  usage  are  also 
computed.  This  also  takes  place  in  subroutine  ACQUIRE.  On  completion 
of  this  task,  control  reverts  to  MAIN  which  then  passes  operations  to 
subroutine  HTLOSS. 

4.  Heat  loss  (Block  6).  At  this  point  the  input  data  file  is  searched 
for  the  maximum  pressure,  the  maximum  time  derivative  of  pressure  and 
selected  values  of  P and  dP/dt.  The  ignition  pressure  is  computed 
using  the  igniter  weight  and  thermochemistry  and  the  theoretical 
(adiabatic)  maximum  pressure  is  calculated.  This  data  is  used  to 

compute  the  Average  Heat  Loss  Rate  (H, ) . Other  operations  including 
setting  the  ignition  time  for  the  analysis,  setting  the  starting  time 
interval,  computing  the  initial  system  temperature  and  gas  generation 
rate  are  also  performed.  These  operations  are  performed  in  subroutine 
HTLOSS.  Control  then  passes  to  MAIN  which  relinquishes  operation  to 
subroutine  SUMUP. 

5.  Run  header  (Block  8).  The  function  of  subroutine  SUMUP  is  to 
provide  the  identifying  information  for  the  run.  Information  on  the 
data  file  designation,  propellant  parameters,  igniter  data,  bomb  data 

as  well  as  summary  information  on  the  pressures  and  the  time  derivatives 
of  pressure  obtained  are  printed.  The  first  page  and  part  of  the 
second  page  of  Appendix  B are  the  output  from  this  operation.  The  only 
portion  of  the  header  not  printed  at  this  point  relates  to  the  initial 
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propellant  surface  area  and  the  number  of  propellant  grains.  Subroutine 
SUMUP  passes  control  back  to  MAIN  which  then  passes  control  to 
subroutine  REDUCE. 

6.  Begin  Computation  (Blocks  9-11).  Subroutine  REDUCE  controls  the 
differential  equation  solver,  subroutine  EVES.  It  specifies  the  number 
of  differential  equations  to  be  solved,  the  stopping  time  and  the 
interval  at  which  data  are  to  be  printed  out.  (See  pages  3 and  4 of 
Appendix  B) . Control  then  passes  from  REDUCE  to  EVES  until  the  stopping 
time  (or  some  other  limiting  condition)  is  reached. 

To  initiate  computations,  subroutine  EVES  calls  subroutine  SETUP 
whose  job  it  is  to  provide  the  initial  values  of  the  functions  and  the 
factors  needed  to  calculate  the  starting  values  of  the  differentials. 
Values  are  assigned  to  the  propellant  mass  generation  rate,  ignition 
aid  burning  rate,  mass  loss  rate,  thermodynamic  parameters,  heat  loss, 
etc.  In  a standard  analysis,  the  initial  heat  loss  rate  is  the  average 
heat  loss  rate  computed  earlier  in  subroutine  HTLOSS.  The  initial  mass 
loss  rate  and  igniter  burning  rate  are  zero.  The  initial  surface  area 
of  the  propellant  is  computed  by  calling  the  appropriate  form  function 
subroutine.  The  surface  area  data  and  the  number  of  propellant  grains 
are  then  printed,  completing  the  "header"  section  of  the  program 
output  (Block  10) . Control  passes  from  SETUP  through  EVES  to  DIFEQ 
where  the  initial  values  of  the  differentials  are  computed.  (The 
functions  evaluated  using  DIFEQ  and  EVES  are  designated  as  Y(N)  and 
their  differentials  as  DY(N);  see  the  FORTRAN  column  in  the  List  of 
Symbols.)  In  a standard  analysis,  the  initial  mass  burning  rate  is 
set  equal  to  the  ignition  aid  mass  burning  rate  at  the  point  of 
ignition.  The  other  mass  differentials  (and  their  integrals)  are, 
of  course,  zero. 

Using  the  starting  values  of  the  function  (Y  ) and  its  differential 
(Y  ) EVES  computes  an  estimai^ed  value  for  the  function  at  the  completion 
or  the  first  time  interval  (Y^ . This  value  is  fed  back  to  DIFEQ  where 

it  is  used  to  obtain  a now  value  of  the  differential  (Yp . The 

difference  between  (Y^)  and  (Yp  is  examined  and  compared  with  an  error 
level  built  into  the  program.  ^If  the  value  passes,  the  value  of  the 
function  (Y.)  and  its  differential  (Y.)  are  accepted,  the  time  interval 
is  incremented  (Block  9}  and  the  process  is  repeated.  The  initial 
integration  steps  are  purposely  toot  small  to  establish  the  initial 
table  of  differences  required  by  the  predictor-corrector  technique. 
Subroutine  EVES  has  the  built  in  capability  to  adjust  the  size  of  the 
time  step  used  as  the  analysis  progresses  so  that  under  conditions  where 
the  predictor  values  are  better  than  required  by  the  difference  criteroon, 
the  time  steps  are  increased  in  size,  providing  a savings  of  computation 
time.  Time  steps  are,  of  course,  also  automatically  reduced  by  the 
program  as  required. 


7.  Compute  Burning  Rates  (Blocks  12  through  22).  The  process  just 
described  is  followed  to  obtain  the  required  values  of  linear  burning 
rates  throughout  the  analysis.  At  each  step  where  a new  value  of  (dw  / 
dt)  is  required,  subroutine  DIFEQ  computes  the  values  based  on  the  ^ 
updated  values  of  the  parameters  in  the  Mass  Burning  Rate  Equation. 

The  value  is  examined  in  EVES  and  accepted  or  rejected  as  necessary. 
Blocks  12  through  22  are  involved  in  the  process;  Of  these  12  through 
19  take  place  in  subroutine  DIFEQ  and  the  appropriate  form  function 
subroutine.  In  a standard  analysis  the  heat  loss  rate  is  constant  and 
the  average  heat  loss  rate  is  used.  Since  propellant  thermodynamic 
characteristics  are  constant,  no  updating  of  values  is  done.  The  same 
is  true  for  both  the  ignition  aid  burning  and  the  mass  loss  terms. 

If  the  solution  of  the  differential  equation  (as  judged  in  EVES)  is 
acceptable  at  any  point,  values  for  burning  rate  are  accepted  along 
with  those  of  w and  dw  /dt.  At  given  intervals  during  the  analysis 
the  data  are  printed  ou?  (EVES  calls  subroutine  PRINT) , see  pages  50 
and  51  of  Appendix  B.  Tests  for  limiting  conditions  (web,  all 
burned,  pressure  = P , and  time  = T ) are  made  throughout  the  analysis 
and  once  one  of  the  limiting  conditions  is  reached,  EVES  returns 
control  to  REDUCE  which  returns  control  to  MAIN. 

8.  Fitting  the  burning  rate  data  (Block  25).  Once  the  reduction 
phase  is  completed,  the  burning  rate  data  between  0. 1-0.8  of  P^  is 
fitted  to  an  equation  of  the  form  roAP".  This  is  done  in  subroutine 
RCALC. 

9.  Plotting  of  data  and  results.  (25,26).  The  question  is  posed  to 
the  operator  whether  the  lower  and  upper  portion  of  the  burning  rate 
curve  are  to  be  suppressed.  Generally  this  is  done,  since  the  most 
meaningful  data  is  between  0.1  6 0.8  of  P . After  this  decision,  P vs.t; 
dP/dt  vs  t;  dP/dt  vs  P/P  and  Inv  vs.  ilnP”aro  prepared  and  the  program 
exits.  The  complete  output  package  from  the  program  may  bo  soon  in 
Appendix  B. 

B.  Nonstandard  Analysis 

If  the  nonstandard  mode  of  analysis  is  chosen,  an  interactive 
display  is  activated  in  subroutine  ACQUIRE  (Block  (4A)  which  allows 
temporary  modification  of  the  pertinent  ballistic  information  appearing 
in  the  data  file.  In  addition  to  those  changes,  decisions  concerning 
a variety  of  options  have  to  bo  made.  A brief  discussion  of  each  of 
the  options,  follows. 

1.  Heat  loss.  (Blocks  7A,  12  and  12A).  In  the  nonstandard  mode,  the 
first  option  concerns  the  assignment  of  heat  loss.  The  standard  option, 
(Block  6).  was  described  earlier.  The  nonstandard  option  (Block  7A) 
is  a much  more  sophisticated  treatment  which  assigns  values  to  the 
convective  and  radiative  components  of  heat  loss  based  on  the  decay 
portion  (after  P ) of  the  firing  record  itself.  The  theory  was 
discussed  earlier.  In  the  reduction  phase  of  the  program 
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the  heat  loss  rate  (H, ) is  computed  at  each  point  in  subroutine  DIFEQ. 

Both  the  convective  and  radiative  heat  loss  rates  are  functions  of  the 
gas  and  wall  temperatures.  Subroutine  WALTEM  is,  therefore,  invoked  to 
provide  current  values  of  wall  temperature  as  are  needed  in  the  analysis. 
This  is  shown  schematically  by  Block  12A  in  Figure  2. 

2.  Propellant  Geometry  (Blocks  18A  S 18B).  The  purpose  of  this  option 
is  to  approximate  the  effect  of  real  propellant  geometry  by  using  a 
distribution  of  web  values  rather  than  nominal  web  values.  The 
computations  take  place  in  DIFEQ  and  the  appropriate  form  function 
subroutine.  A web  deviation  value  may  be  input  in  the  nonstandard  mode 
resulting  in  the  proportioning  of  the  propellant  charge  into  five  parts: 
two  each  of  10  percent,  two  each  of  20  percent,  and  one  each  of  40  percent 
of  the  total  charge  with  webs  differing  from  the  nominal  by  plus  or 
minus  suitable  factors  times  the  web  deviation  derived  from  a normal 
population  distribution  curve.  These  five  propellants,  equivalent 

in  weight  to  but  not  having  the  same  initial  area  or  the  same  overall 
form  function  as  the  total  charge  considered  as  a nominal  case  are 
treated  as  separate  entities  in  the  burning  area  portion  of  the 
reduction  and  the  instantaneous  areas  summed  for  calculation  of  the 
linear  burning  rate.  This  option  eliminates  sharp  discontinuities  of 
surface  area  at  certain  points.  (Computer  generated  pressure-time 
curves  using  this  approacli  have  resulted  in  much  closer  agreement  with 
dP/dt  data  taken  from  real  propellant  geometries.)  Only  the  burning 
area  portion  of  the  analysis  is  affected. 

3.  Ignition  Deviation.  This  option  permits  the  simulation  of  flame 
spreading  effects  that  may  occur  in  on  experimental  firing.  A 
population  distribution  similar  to  the  one  above  (propellant  geometry 
section)  is  used.  The  propellant  is  assumed  to  be  composed  of  five 
samples,  each  ignited  at  a different  time.  (It  is  the  value  of  the 
time  interval  that  is  input  from  the  keyboard.)  Hie  distance  burned 
for  each  fraction  is  calculated  separately.  Those  values  of  x are 
used  to  compute  the  area  of  the  burning  surface  of  the  d>arge  during 
the  analysis.  A linear  interpolation  method  is  used  to  ’'emove  gross 
discontinuities  In  the  burning  area  between  the  times  of  ignition  of 
two  adjacent  samples  and  upon  burnout  of  the  first  sample.  Only  the 
burning  area  portion  of  the  program  is  affected  by  this  option. 

Operations  take  place  in  subroutine  DIFEQ  which  calls  the  appropriate 
form  function  subroutine  as  required. 

4.  Ignition  Aid.  (Block  IS,  ISA),  ITus  option  is  useful  in  describing 
situations  in  which  not  only  an  igniter  but  also  an  ignition 

aid  is  present.  Since  the  ignition  aid  actually  has  a finite  bum  time 
relative  to  the  propellant,  the  option  allows  modification  of  the 
"ignition  pressure"  to  account  for  the  partial  burnout  of  the  aid 
material  at  the  tine  that  propellant  ignition  takes  place.  During  the 
early  pJiase  of  the  firing,  therefore,  both  propellant  and  ignition  aid 
will  bo  contributing  to  the  mass  generation  rate,  hence  to  dP/dt.  The 
aid  contribution  ciust  bo  evaluated  separately.  At  present  this  is 


accomplished  by  external  input  or  by  analysis  of  the  dP/dt  portion  of 
the  record  just  prior  to  the  chosen  ignition  point  (pressure),  translat- 
ing this  dP/dt  into  an  equivalent  mass  burning  rate  of  the  aid  and  fitting 
it  to  a power  law  curve.  Evaluation  of  the  instantaneous  value  of  the 
ignition  aid  burning  rate  is  performed  in  subroutine  DIFEQ.  The 
function  is  integrated  in  EVES.  The  computation  of  the  mass  burning 
rate  (Block  17)  is,  of  course,  affected  (see  Equation  14  term  E) . 

5.  Variaoie  thermochemistry.  (Blocks  14  and  14A) . This  option  allows 
the  input  of  the  required  thermochemical  information  for  the  propellant 
in  the  form  of  a table  of  values  versus  pressure.  The  table  entry 
takes  place  in  subroutine  ACQUIRE.  Values  for  the  rate  of  change  of 
covolume,  heat  capacity  etc.  are  computed  in  subroutine  DIFEQ  and 

used  in  calculating  the  mass  burning  rate  of  the  propellant  (also  in 
DIFEQ) . 

6.  Vented  Chamber  Operation.  (Blocks  16  and  16A) , This  option  allows 
analysis  of  firings  made  with  '“leaking''  or  vonted  vessels.  The 
effective  vent  area  and  the  blow-out  pressure  of  the  vent  are  necessary 
inputs  and  their  presence  will  automatically  result  in  analysis  beyond 
the  time  of  P^^^  and  the  computation  of  the  mass -discharge  term  in  the 

differential  equations.  Instantaneous  temperatures  as  well  as  the 
effects  of  changing  molecular  weight  are  used  in  determining  this  term. 

7.  Burning  Surface.  (Block  ISA).  In  this  option,  an  arbitrary  total 
burning  surface  area  as  a function  of  some  characteristic  burning 
dimension  nay  bo  input  to  the  ?i*ogram  if  none  of  the  available  foim 
functions  appear  suitable.  The  input  consists  of  a table  of  points, 
and  the  instantaneous  area  is  interpolated  linearly  from  the  information 
furnished.  For  seldom  used  unique  geometries,  it  avoids  the  nuisance 

of  rccom|>iling  the  program  to  include  a special  form  function  generator. 

C.  Summary 

A versatile  closed  bomb  data  reduction  program  has  been  developed 
to  compute  linear  burning  rate  of  propellants  from  the  pressure-time 
histories  of  their  burning  process.  In  .Section  II,  the  capability  of 
the  program  ss  seen  from  the  user's  point  of  view  was  discussed,  and  a 
comparison  of  the  standard  and  nonstandard  modes  of  analysis  was  made. 

An  introduction  to  the  theoretical  treatment  was  made  in  Section  III 
and  the  relationship  between  the  equations  used  in  the  program  and  the 
simplified  derived  equation  was  examined.  In  Section  IV,  program 
structure  and  operation  were  ©.xamined  and  a flow  chart  of  the  program 
as  well  as  a diagram  of  the  subroutine  hierarachy  was  given.  Finally, 
the  impact  of  each  of  the  options  on  program  execution  was  discussed. 

The  derivation  of  the  mass  burning  equations  (Appendix  0)  as  well  as 
a program  listing  (Appendix  G)  are  provided. 
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APPENDIX  A 


Input  Requirements  for  CBRED 


Input  Requirements  for  CBRED* 


1.  Run  ID.  Provides  input  file  designation  for  program. 

2.  Propellant  Data. 

a.  Type 

b.  Weight  (grams) 

c.  Density  (grams/cc) 

d.  Grain  Type 

e.  Length,  OD,  ID  (inches) 

f.  Inner  Web,  Outer  Web  (inches) 

g.  Theoretical  Impetus  (ft-lb/lb) 

h.  Flame  Temperature  (•K) 

i.  Average  Molecular  Weight  of  Products  (grams /mole) 

j.  Covolume  (cubic  inch  per  pound) 

k.  Gamma  (ratio  of  specific  heat) 

3.  Igniter  Data 

a.  Type 

b.  Weight  (grams) 

c.  Impetus  (ft-lb/lb) 

4 . Experiment  Data 

a.  Bomb  Volume  (cc) 

b.  Bomb  Temperature  (•’’,) 

* The  input  information  lifted  is  actually  what  is  required  by  the 
program  being  documented.  Conversion  to  metric  units  is  being 
planned. 
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5.  Data  Arrays  (digital) 

a.  Pressure-Time*  (Kpsi,  milliseconds) 

b.  Time  Derivative  of  Pressure-Time*  (mega  psi,  milliseconds) 


*The  input  information  listed  is  actually  what  is  required  by  the 
program  being  documented.  Conversion  to  metric  units  is  being 
planned. 
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APPENDIX  B 


Sample  Output,  CBRED* 


•Program  output  is  currently  not  in  metric  units.  Conversion  is  being 
planned. 
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Eqtiation  for  FORM  SP 


Form  Function  of  a Sphere. 
The  Distance  Burned 


S,  Surface  Area,  as  a function  of  x, 


W * D 

where:  D = initial  diameter  of  sphere 
W a propellant  web 
X = depth  burned  at  time  t. 

S = 0 for  W ^ 2x 
Otherwise : 

D'  = D-2x 
Surface  Area 


Equations  for  FORM  CY 


Form  Function  of  a Right  Circular  Cylinder.  S,  Surface  Area  as  a 
Fvmction  of  x,  the  Distance  Burned 

W = D 

where:  D = initial  grain  diameter. 

W = propellant  web 

X = depth  burned  at  time  t 

S = 0 for  L _<  2 X 

(L  = Initial  grain  length) 

S = 0 for  W < 2 X 

Otherwise: 

D'  = D-2x 
L'  = L-2x 


End  Area: 

E « Tr/4  (D*)^ 

Surface  Area: 

S = 2 E+ir  L’  D' 
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Equations  for  FORMl 


Form  Function  of  a Single  Perforated  Right  Circular  Cylinder  (Axially 
Symmetrical).  S,  Surface  Area,  as  a function  of  x,  the  Distance 
Burned. 


where : 

W » propellant  web 


D = initial  grain  diameter 

Dp  = initial  perforation  diameter 

X * depth  burned  at  time  t. 

3 = 0 for  L 2x 

(L=  initial  grain  length) 

3=0  for  W < 2x 


Otherwise : 


D'  = D-2x 
L'  = L-2x 


Dp  ♦ 2x 


End  Area: 

E =1  [ (D’)2  - (0  )2) 
4 P 


Surface  Area; 

S=2E  ♦ itL*  (O’  ♦ Op,) 
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Equations  for  FORM  7 


Form  Function  of  a Seven  Perforated  Right  Circular  Cylinder  (Axially  . 
Symmetrical).  S,  Surface  Area,  as  a Function  of  x,  the  distance  burned. 

I,  To  splintering; 


W = D - 3 D 

E 

4 


where : 


W = propellant  web 

D = initial  grain  diameter 

Dp  = initial  perforation  diameter 

X = depth  burned  at  time  t 

S = 0 for  L'  £ 2x 

(L'  = instantaneous  grain  length) 

S = S for  W=2x 
max 

S = < S for  W<2x 
max 

for  0 £ X,  define 

D'  = D-2x 
L'  » L-2x 

(L  = initial  grain  length) 

D , =*  D ♦2x 

p‘  p 

Then,  for  0^2x^W 
End  Area 

E » 1 [ (D*)^  - 7 (D  ,)2  3 
4 P 

Surface  Area: 


S = 2E  ♦ It  L'  (D'  ♦ 7 Dp,) 
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Form  Function  of  a Seven  Perforated  Right  Circular  Cylinder  (Axially 
Symmetrical) . Surface  Area,  S,  as  a Function  of  x,  The  Distance  Burned 
(continued) . 


II.  After  Splintering* 

define:  W.  »D„+W  and  let 
w p 

C.min  ll.  D^-D 

( 2(D.Dp-W^^)  I 

Then  for  W<2x<C,  let 


= End  area  of  outer  slivers 
for  ct  < TT  /6 

Ej  « 3 D'  W^sin  a + 3/2  J^(U')^(n/6-o)-W^^  ^ -(Dp,)^(Vl/2  sin  Q)j 

Ej  = 0 for  a^»/6 

° Surface  area  of  outer  slivers 
for  a<  tr/6 

= 2Ej^(60Dp'*(ff-6a)  D')  l' 

» 0 for  a /6 

* Treatment  developed  by  Hr.  Franz  Lynn,  USABRL. 


6S 


E2  = End  area  of  inner  slivers 
for  0<tr/3 


3/2 


[■-’  f - 


3/2 (Dp,)' 


(sin  0+  tt/3  - 0) 


- 


E2  = 0 for  02^ir  /3 

$2  = Surface  area  of  inner  slivers 
for  0 < TT  /3 


S2  = 2E2  + 9Dp,  (ir/3-0)L' 

S2  “ 0 for  TT  / 3£0 

S = Total  surface  area  of  slivers 
S = + S2 

E = Total  end  area 


for  C < 2 X,  EaO,  S=0 


Equations  for  FORM  19 


Form  Function  for  a Nineteen  Perforated  Right  Circular  Cylinder 
(Axially  Symmetrical) . S>  Surface  Area,  as  a Function  of  x.  The 
Distance  Burned. 

I . To  Splintering 

W = D-5  D 

E 

6 


where : 


W = propellant  web 

D = initial  grain  diameter 

D = initial  perforation  diameter 


X » depth  burned  at  time  t 

S = 0 for  L'  < 2x 

(L'  = instantaneous  grain  length) 

S = S^„^  for  W • 2x 

S < s"*®*  for  W<2x 
max 


for  0;<  X,  define: 


D'  = D-2x 
L»  = L-2x 

(L»initial  grain  length) 
Dp^  = Dp  + 2x 

Then,  for  0 _<2x^W 


End  Area 


E l(D*)^  - 19  (D  .)^  ] 
4 P 


Surface  Area 


S 2E  L*  (D*  + 


19  Dp.) 
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Fc  ' Function  of  a Nineteen  Perforated  Right  Circular  Cylinder  (Axially 
Symii,  trical) . Surface  Area  S,  as  a function  of  x,  The  Distance 
Burned,  (continued) . 

II»  After  Slivering* 

.*■  • 

define:  W = + W and  let 

w p 


C = min  /L 


12[(D^Dp)^  - 16W^^] 
[D+Dp]2  - 12W^2 


Then,  for  w<2x<C,  let 


0=2  cos'^  ) min  , 1 


aj  = cos’"  I min((l/8  [(D')^  - ♦ 2 


<^2  = cos 


-1 


= cos 


-1 


I min((l/4  [(D')^].  (Dp,)2  - (D^,)^]  * 

I max((l/8  [(Dp,)2  - (D')^]  + 2 W„^)  ^ 

Vp-  ’ /) 


&2  ° cos 


-1 


I max^( 


(1/4  [(Dp.)^  - (D')’]  >5  W„^) 


“A'  f 


and; 


a=ai  ^02 

B=Bl  *&2  • 0 " STr/6 

• Treatment  developed  by  Mr,  Franz  Lynn,  USABRL. 
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yp 

if*’- 

Jv- 


a-' 

.4^,. 


= End  area  of  Outer  Slivers 
for  o<  IT  /6 


E.  = 3 D'  W (2  sin  a,  ♦ 
1 w 1 


^ sin  02)  - 6 


+ 3/2  [(D')^  (it/6-tt)  - (Dp,)^  (sin  9 + 0)] 


Ej^  = 0 for  a^ir  /6 


= Surface  Area  of  Outer  Slivers 
for  a<nr  /6 

S^  = 2E^  .+  6 [ D*  (ir/6-o)  + D^,  B]  • L* 

Sj^  = 0 for  a^'ir/6 

E2  = End  area  of  Inner  Slivers 
for  0<  IT  /3 

E,  = 6 [wj  JT  - 3/2  (sin  0in /3  - 0)} 

z w ^ p' 

£2  = 0 (or  0>  i»  /3 

S^  = Surface  Area  of  Inner  Slivers 
for  0<  » /3 

S,  = 2E,  ♦ 36  D ,(if/3  - 0)L' 

— 4 p 

S^  = 0 for  0^w  /3 

S = Total  Surface  Area  of  Silver's 
S S^  . S2 

E a Total  End  Area 

C . U;  . Cj 

for  C<2x,  U = 0 and  S = 0 
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APPENDIX  D 

Derj.vation  of  Mass  Burning  Rate 
Equation  for  CBRED 


Derivation  of  the  Mass  Bitrning  Rate 
Equation  for  CBRED 


The  same  order  of  presentation  is  followed  as  in  the  text.  The 
Equation  of  State  is  presented  first,  the  Energy  Equation,  second, 
followed  by  the  Mass  Burning  Rate  Equation. 


(1)  Equation  of  State  of  Gas 


PV  = w R 
s s s 


T 

s 


(23) 


where : 


P = pressure 
V = system  volume 


(c. 


c ) 
_£ 


w ) 


- w n 
s ' 


= empty  bomb  volume 
c a starting  weight  of  ignition  aid 
c a starting  weight  of  propellant 

V a solid  propellant  density  (assumed  same  for  ignition  aid) 
“g  a weight  of  ignition  aid  burned* 
w a weight  of  propellant  burned* 

"s  "r.VitV^V  , 
w^  = weight  of  air  in  chamber 

w*  a weight  of  initiator  combustion  products  in  chamber 

ignition  aid  combustion  products  in  chamber* 
w^S  weight  of  propellant  combustion  prouucts  in  chaabor* 
covolume 


R R 
s 


u 

W . 


Ry  * uiviversal  gas  constant 


o 


ai 
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Mj.  a molecular  weight  of  air  (taken  as  29.) 

Mi  a molecular  weight  of  initiator  combustion  products 
Mg  = molecular  weight  of  ignition  aid  combustion  products 
Mp  a molecular  weight  of  propellant  combustion  products 

* Note:  For  a closed  bomb  system  there  is  an  </bvious  redmdancy  between 
Kp  and  Wpi  and  Wa  and  Wgi.  But  for  a leaking,  or  vented  vessel, 
tne  distinctions  are  important. 


(2)  Energy  Balance  Equation 

It  is  more  convenient  to  describe  the  energy  balance  dynamically. 
The  following  governing  equation  applies: 


t T„  w ] - H,  -C.  T 
a Op  p-*  IPs 


dt 


(24) 


where : 


Cy  - heat  capacity  at  constant  volume  (assumed  same  for 
ignition  aid) 

Tq^=  isochoric  adiabatic  flame  temperature  of  the  ignition  aid 
Tqp=  isochoric  adiabatic  flame  temperature  of  the  propellant 

w = AP^  By  definition.  Ignition  aid  mass  burning  rate. 

d. 


V’ 


p 


= mass  burning  rate  of  the  propellant. 


= beat  loss  rate 

Cp  = heat  capacity  at  constant  pressure 
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where : 


"n  ” * St  \ 

£,  * gravitational  constant 
Pst  “ stagnation  pressure 

= effective  throat  area  (sonic  control  assumed) 
C*  = characteristic  discharge  velocity 


Tst  ® stagnation  temperature 
Ms  system  molecular  weight 
2 

^ =0  function  of  the  specific  heat  ratio 


75 


(3)  Rate  of  Conversion  of  Solid  to  Gas 


Solving  the  Equation  of  State  (23)  for  and  differentiating  yields: 


dT^  1 
s 

dt 

where ; 


R w 
s s 


PV  +V  P - PV  (w  R + R 
s s s s s s 

. R w 

s s 


w ) 

. s , 


(25) 


V^=-w^  n-,ni^  (w^) 


n = ^ ^ By  definition 

dP  ^ dt 


R 


nij,  - nip 


W 


ni_  = w w.  w,  w w,w 

T _£.■*■_!  * P - pi  n 

M M.  M M M w 


I a 


P s 


W 

-£LE 

(Mp)' 


w 8 - w X w By  definition 
r r n 

W 


w,  a - w,  X w By  definition 
i n 

w 


w , o w„  - w , w By  definition 
al  a ajL  n '' 
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w 


w + 
r 


''i  '*al 


* w_ 


_El 

w 


Note  also: 


w - w , w By  definition 
P ^ 


Differentiating  the  left  hand  side  of  Equation  (24)  as  indicated  gives; 


C.,  w T + C.,  T w + w T C,,  » C,,/T»  w + T.  w 
vss  vss  ssv  vl  Oa  a Op  p 


(26) 


-H,-  T w 
L p s n 


where : 


Cy  - d (Cy) 
dP" 


dP  by  definition 

at 


Solving  Equation  (26)  for  the  rate  of  change  of  system  temperature 
with  time  (T^,)  one  may  simultaneously  solve  with  Equation  (25) , which 
on  expansion^of  terms  and  appropriate  manipulation  yields  the  following 
equation  for  mass  burning  rate: 


O 
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APPENDIX  E 


Equations  for  Wall  Temperature 
Computations  in  WALTEM 


Heat  loss  is  computed  by  coupling  transient  conduction  in  the  bomb 
wall  with  the  convective  and  radiative  transfer  to  the  wall.  In  the 
wall,  the  governing  equation  is 


6T  5^T 
w w ^ 


(23) 


where  « temperature  field  within  the  wall 

a » thermal  diffusivity 

X = radial  distance  into  the  wall 

The  initial  condition  is  a uniform  temperature  (set  at  298K) . The 
boundary  conditions  are: 

atX.O  : 

at  X . . |I  = 0 


An  explicit  centered  difference  scheme  is  used  for  the  calculation.  At 
the  boundaries  the  scheme  is  as  follows: 

at  X = 0 : 


q&t 

(ax)‘ 


2 (AX) 


- 2Tj  ♦ 2T 


n 

2 

Ml 


whore  n ♦ 1 is  the  now  time  level 

n the  old  time  level 

Tj  the  first  interior  point 

^2  the  second  interior  point 


at  back  boundary  X ■ N AX: 

O 


T = T , 

U H -1 

N * number  of  grid  points 
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appendix  f 


Capsule  Su»ary  of  Suboutlnes 


S3 


Capsule  Sunaary  of  Subroutines 


ACQUIRE 

HTLOSS 

WTFIT 

SUMUP 

REDUCE 

EVES 

SETUP 

DIFEQ 

PRINT 

F INDTP 

FINOl 

SIHEQ 

RCALC 

VPLOT 

Registered 


Subroutine  to  get  file  data,  allow  update  and 
start  process 

Subroutine  to  calculate  average  heat  loss  value, 
average  heat  transfer  coefficient 

Subroutine  to  fit  the  decay  portion  of  the  P-t 
curve  to  find  heat  loss  coefficients. 

Subroutine  to  print  summary  sheet  of  analysis. 

Driver  for  the  differential  equation  solver. 

Subroutine  solves  N simultaneous  first  order 
differential  equations  by  the  Adams  method. 

Subroutine  sets  the  initial  conditions  for  the 
integration  and  defines  the  accuracy  required  in 
the  solution. 

This  subroutine  evaluates  the  derivative  values  at 
each  proposed  step  in  the  solution.  On  each  call 
Tl  contains  the  current  time,  and  Y's  are  the  current 
integrated  values.  The  results  of  a call  may  be 
rejected  and  the  step  size  cut  in  order  to  maintain 
accuracy . 

This  subroutine  is  called  to  output  accepted 
values  during  the  integration  procedure. 

A subroutine  for  table  lookup  with  linear  interpolation 
A direct  access  read  lookup  modification  of  FINDl. 

A subroutine  for  tabic  lookup  with  linear  interpolation 
Extrapolates  values  out  of  tabic  range,  remembers 
last  argument  value  and  begins  search  from  that 
value. 

A subroutine  solving  L equations  in  M unknowns,  with 
N sots  of  right  hand  constants. 

Subrout^'fio  to  do  least  squares  fit  on  bun^  rate  data. 
Vorsatoch*  uiiit  plotting  package. 


trademark,  does  not  constitute  endorsement  by  the  US  Army. 
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APPENDIX  G 
* 

Program  Listing . C6RGD 


RT-U 


0601 


0602 

0603 

6604 

0605 

0606 
0637 
0600 
6609 
0010 
0011 
0312 

0014 

0015 
0616 

0017 

0018 
8019 
6020 
0021 
0022 
0023 
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C MAIN  DRIVER  FOR  LATEST  REDUCTION  PROGRAM.  10/14/75 

COmON  RUNID<5) .RTITLEt IB) .DATE(3) .0PERN(5)  .PROPNC 10) . 

1 PSORCEC 10) .PRLOT( 10) ,PREM( 10) .TIGNRC 10) .GAGE( 10) . 

2 SPACE1(36).CHGIJT.L1ID,PTEMP,UGHTI.PCORR.BVOL,0TEMP. 

3 SCfiL.CALCQN.STIME.FFID.GL.OD.PD,UOD,Fl.DUM.XMJXl.ETl.GAMl. 

4 SPACE2(4). ISKl.KTOT.MZ.MY.PMAX.TMAX.NPMA.PPM.DPMAX.  IHL.H, 

5 RHOC.RHO.T1Q.T90.T1090.P9(5).DP9(5).NPTH.PX(20). 

6 FX(20).ETX(20).XMUX(20).GAMX(20).IDFF.NPA.IJEBX(20).ABX(20). 

7 ISK, ISK2,MEti2.XMl.RP.TZERD,V0L.CP,CST,TM,PIGH.TIGN.PTHE0. 
a C0Nl.CQN9.MEM3.MEM7,SFA2.L9.P,DP.TST0P. 19, 

9 LB. CONU, TBU, FRAC. 2.SRAT. XTBU. RATE. RL 
COmiN  /GENE/  UIDP(5).LJ0DP(5).UCR(5).SFAC<5.S).TDC!.(5). 

I UDEV.TDEV.TUST 

COMMON  /BLAH/  PLO,  PHI.  ACO.  XNC.  RSQ.  RMS 
C0M10N  /KEVIN/  TE(35).  XT(35) 

DATA  DUm/'UQU  V 
URITE  (11.197)  UQU 
197  FORMAT  (A4) 

REUIND  11 
CALL  ACQUIRE 
CALL  HTLOSS 
SPACE2(3)  - 0.8 
IF  (IHL.EQ.2)  GO  TO  1976 
CALL  HTFIT 
1976  CONTINUE 
CALL  SUMUP 
CALL  REDUCE 
CALL  RCALC 
TtIPE  100 

100  FORMAT  (///) 

CALL  VPLOT 

STOP 

END 
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8801  SUBROUTINE  ACQUIRE 

C SUBROUTINE  TO  GET  FILE  DATA.  ALLOU  UPDATE  AND  START  PROCESS 
0002  COMMON  RUNID<5).RTITLE(18),DATE(3),0PERNC5).PR0PN(18). 

1 PSORCE C 10) ,PRL0T< 18) .PREM( 10) , TIGNR ( 10) .GAGE ( 10) . 

2 SPACE  1(36). CHGUT. UID. PTEhP. UGHTI . PCORR.BVOL,  BTEMP. 

3 SCAL.CALCON.STIME.FFID.GL.OD.PD.UOD.Fl.DUM.XMUXl.ETl.GAm. 

4 SPACE2(4).ISKl.KT0T.ri2,t1Y.PMAX.TMAX.NPW.PPM.DPMAX.  IHL.H. 

5 RHOC.RHO.T10.T9a.T109O.P9(5).DP9(5),NPTH.PX(20). 

6 FX(20).ETX(20).XrtJX(20).GAMX(20).IDFF.NPA.UEBX(20).RBX(20). 

7 ISK. ISK2,MEM2.XMI.RP.T2ERO.VOL.CP.CST.TM.PIGH.TIGN.PTHEO. 

8 CONI. C0N9.MEM3. MEM?. SFAZ.L9.P.DP.TST0P. 19. 

9 L8.C0NU.TBU.fr AC. 2. SR AT. XTBU 


0003 

COltWN  /GENE/-  UIDP<5).UQDP(5).UCR(5).SFAC(5.5).TDEL(5) 
1 UDEV.TDEV.TUST 

0004 

COMiON  /BLAH/  PLQ.  PHI.  ACO.  XNC 

0005 

DIf-ENSION  A9(4) 

0006 

DATA  AY/'Y'AA9/0.25.  .375..S0..62S/'.BL/'  'A 

l.AVGAAVE  'API/' IP  'AP7A7P  'AP19A19P  'A 

2 SP/"SPH  'ACYACY 

0007 

DIIENSION  LIST  (400) 

8000 

EiSUIVALENCE  (L IST(  1) .RUNID(  1) ) 

0009 

EQUIVALENCE  < IBUG.LIST(390) ) 

0010 

EQUIVALENCE  (ISK1.LIST(399)>. (KTOT.LIST(4B0)) 

0011 

EQUIVALENCE  ( ATH.  SPACE  1(10)).  (PBLOU.  SPACE  KID), 

1 (FA  ID. SPACE  1(12)). (TAID.SPACEI ( 13) ) . (CAID. SPACE  1(15)) 

2 (XMA.SPACEKU)) 

0012 

CALL  ASSIGN  (2. 'SYsCBDAT.OOl'.O. 'SCR') 

0013 

DEFINE  file  2 (3O0Q.6.U.M2) 

0014 

TYPE  300 

0015 

300 

FORINT  (SX. 'ENTER  SMOOTHED  TAPE  ID'/) 

8916 

CALL  ASSIGN  (3. 'FILNAM.EXT'.-l.'RDO'') 

8017 

DEFINE  FILE  3(3l00.4.U.t1Y) 

00  le 

UDEV  • 0.0 

0819 

TDEV  • 0.0 

0020 

DO  69  1 • 1.130 

0021 

J - 4*1 

0022 

L « J - 3 

0023 

69 

READ  (3*1)  (LIST(K).  K * L.J) 

0024 

K « 0 

0025 

IfiUG  - 0 

8026 

DPMAX  • 0.0 

002? 

RHOC  - SPACE2(l) 

0020 

PMAX  - 0.0 

0029 

PAID  - 362080, 

0050 

TAID  - 3080. 

0031 

RU  - lS44.»»<1.0 

0032 

KMA  •>  RU*TAID/FAID 

0033 

CAID  • 0.0 

0834 

ATM  - 0.0 

0333 

PBLOU  -0.0 

3936 

KM  - KTOT  + 180 

0037 

DO  68  J • 101. KM 

90 


0038 

READ  (3'J)  01,  Cl 

0039 

I - J - 100 

0040 

11  - <l-l)i«<ISKl  + 

0041 

TI  - n’«TIME>t<l.0E 

0042 

K ■ K + i 
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0843 

IF  (8l.LT.PmX)  GO  TO  88 

0045 

pmx  - Bl 

0046 

Tmx  - TI 

0047 

NPm  - K 

6048 

86 

IF  (Cl.LT.DPmX)  GO  TO  68 

0050 

PPM  ■ Bl 

005  i 

DPmX  - Cl 

0052 

68 

URITE  (2*1)  TL  Bl,  Cl 

0053 

MEM?  - 1 

0054 

PPES  - 0.99it^mX 

8055 

CALL  FiNDTP  (PBES,TMAX,3>2.  l.2,NPMfi,(«M7,KZ) 

0036 

PDES  - a.wmx 

0057 

CALL  FIKOTP  <P!)ES,T10,3,2,  l,2,NPm,MEM7,r«) 

8058 

PD£S  - 0.3>i<PMAX 

0859 

CALL  FINDTP  (PBES,T90,3,2.  l»2»NPm,l€M7,MZ) 

0860 

ri098  - T98  - Tie 

08SI 

SUM  8.8 

9062 

DO  62  I - 1,4 

0063 

PSd)  - A9(l)*PmX 

0064 

CALL  FINDTP  CP3(I),DP9<l),3,2,3.2,NPm,f1EM7, 

0865 

62 

SUM  • SUM  4-  DPS  (I) 

6066 

DPS (5)  • SUrV4.0 

0867 

P9(5)  - AVG 

0068 

ENDFtLE  3 

0069 

FRAC  « 1.0 

0073 

IHL  • 2 

00?' 

H - 0.0 

0072 

XMI  - 60.62 

0073 

IF  (PCORR.GT. 150880.)  XMI  - 25. 

8075 

TYPE  381 

0076 

381 

FOWtAT  (SX.'IS  THIS  TC  BE  A'v' 

1 SX. 'STANDARD  ANALYSIS,  Y OR  N7*) 

897? 

ACCEPT  302.  AHSP 

0070 

382 

FORMAT  (Al) 

8079 

CALL  PLOT  (27.12,4) 

8830 

SPACE2(2)  - 10.1 

8501 

FLO  - 0.1 

8682 

PHI  « 9.9 

8893 

IF  (ANSP.EQ.AY)  GO  TO  39 

6885 

TYPE  303,  RHOC 

8086 

303 

FOSmT  (5X. 'DENSITY*, F12.5) 

0607 

ACCEPT  364.  TEff 

8086 

304 

FORMAT  (66 12.0) 

8889 

IF  (TEm.HE.e.O)  RHOC  - TEW 

8091 

TYPE  385.  PCORR 

8892 

385 

FORMAT  (SX.* IGNITER  IMPETUS*. F12.S) 

6893 

ACCEPT  304.  TEMP 

6034 

IF  (TEMP.NE.B.8)  PCORR  - TEM* 

0898 

TYPE  366.  UCHTl 

0897 

format  (SX.* igniter  UEIGHT*.F12.S) 

8036 

ACCEPT  304.  TEm 

0839 

IF  <TEm.HH.0.0)  UGHTI  - TEMP 

92 
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0iea 

TYPE  1305,  PAID 

0109 

1305  FORMAT  (SX,* IGNITER  AID  IMPETUS' >F 12. 5) 

0110 

ACCEPT  304.  TElf» 

0111 

IF  (TEMP, NE. 0.0)  PAID  - TEMP 

0113 

TYPE  1306.  CAID 

0114 

1306 

FORMAT  (5X,' IGNITER  AID  UEIGHT'.F12.5) 

0115 

ACCEPT  304.  TEMP 

0110 

IF  (TEMP. NE. 0.0)  CAID  ■ TEMP 

0118 

CAID  - CAID.'454. 

0119 

TVFE  1396.  TAID 

0120 

1396 

FORMAT  (5X.' IGNITER  AID  TEMPERATURE' .F 12. 5) 

0121 

ACCEPT  304.  TEMP 

0122 

IF  (TEMP. NE. 0.0)  TAID  - TEM» 

0124 

XMA  - RU>KTAID/FArD 

0125 

CALL  PLOT  (27.12.4) 

0126 

TYPE  307.  FFID 

0127 

307 

FORMAT  (5X. 'GRAIN  TYPE'.2X.A4) 

0128 

ACCEPT  308.  TEhf> 

0129 

308 

FORMAT  (A4) 

0130 

IF  (TGMP.NE.BL)  FFID  - TE^P 

0132 

TYPE  339.  H 

0133 

309 

FORMAT  (5X,'HKAT  LOSS  NUMBER' .F12. 5) 

0134 

ACCEPT  304.  TEMP 

0135 

IF  CTEMP.NE.0.0)  H - TEMP 

0137 

TYPE  310.  IHL 

0138 

310 

FORMAT  i5X.'H£AT  LOSS  OPTION'. 14) 

0139 

ACCEPT  311.  IT 

0140 

311 

FORMAT  (ID 

0141 

IF  (IT.NE.0)  IHL  - IT 

0143 

TYPE  991.  SPACE2(2) 

0144 

91  i 

FORMAT  (5X,'B0MB  UALL  AREA'.  F12.5) 

0145 

ACCEPT  304.  TEMP 

0146 

IF  (TEMP. NE, 0.0)  SPAC.E2(2)  • TEMP 

0148 

TYPE  312 

0149 

312 

FORM/vr  (5X. 'AVERAGE  THERMOCHEMISTRY'/ 

1 SX,  IS  TO  USED.  Y OR  N7') 

0150 

ACCEPT  302.  ANST 

0151 

CALL  PLOT  (2M2.4) 

0152 

IF  (ANST.EQ.AY)  GO  TO  90 

0154 

TYPE  313 

0155 

313 

FORMAT  (SX.‘ ENTER  NUMBER  OF  TABLE  ENTRIES') 

0155 

ACCEPT  311.  NPTH 

3157 

lYPE  314.  NPTH 

0158 

314 

FORMAT  (SX. 'ENTER  PRESSURE,  IMPETUS. COVOLUME'/ 

1 5X.' FLAME  TEMP  AND  SPECIFIC  HEAT'.2X, I 1.2X,' ENTRIES') 

0159 

DO  44  I « l.MPTH 

0160 

44 

ACCEPT  304.  PX(I).  FX(I).  ETX(D.  XMJX(I),  GAMX(I) 

0161 

RUl  - 1544.)P1.8 

0162 

RU2  - RU 1/778. 

0163 

XMUXl  - RUI*XMUX(1)/FX<1) 

0164 

GAMl  - RIJ2/(XMJXl*GAMX(l))  + 1.0 

0165 

GO  TO  97 

94 


01GS  98  TYPE  323,  FI 

0167  323  FORMAT  (5X, 'PROPELLANT  IMPETUS', F 12. 5) 

0160  ACCEPT  304,  TEMP 

0169  IF  (TEMP. NE. 0.0)  FI  • TEhP 

0171  TYPE  324,  ETl 
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0t?2 

0173 

0174 

0175 
0177 
0170 

0179 

0180 
0181 

0183 

0184 

0185 

0186 
0188 

0189 

0190 

0192 

0193 
0195 
019? 

0198 

0199 

0200 
0201 
0202 

0203 

0204 

0205 

0206 

0207 

0208 

0209 

0210 
0211 
0212 
0214 
0216 
0218 
0220 
0222 

0224 

0225 

0226 


0227 

0228 

0229 

0230 


324  FORMAT  (SX, 'PROPELLANT  COVOLUME' <.F12. 5) 

ACCEPT  304,  TEMP 

IF  (TEMP. ME. 0.0)  ETl  - TEMP 
TEMPI  - XMUXl 
TEMP2  - GAMl 
TYPE  325,  XMUXl 

325  FORMAT  (5X, 'FLAME  TEMP ',F 12. 5) 

ACCEPT  304,  TEMP 

IF  (TEMP. ME. 0.0)  XMUXl  - TEMP 
TYPE  326,  GAMl 

326  FORMAT  (5X.' SPECIFIC  HEAT*, FI 2. 5) 

ACCEPT  304,  TEMP 

IF  (TEMP, ME. 0.0)  GAMl  - TEMP 
TEMPI  - XMUXl 
TEMP2  - GAMl 

IF  (XMUX1.lt.  1000.)  TEhPl  - F1>kXMUX1/(  1344.*! .8) 

TEMP3  - XMUXl 

. IF  (XMUXl. GT.  1000.)  TEMP3  - 1544.*1 .0»XMUX1,'P1 

IF  (GAMl. GE. 1.0)  TEMP2  - 1544. ♦I .8/778. /TEMP3/(GAM1  - 1 
NPTH  - 4 
P2  - 1.0 
DO  447  1-1,4 
PX(I)  - P2 
FX(I)  - FI 
ETX(I)  - ETl 
XMUX(I)  • TEMPI 
GAMX(I)  • TEMPO 
P2  - P2  + 1.0 
447  CONTINUE 

RUl  • 1544.'e<1.8 
RU2  - RU 1/778. 

XMUXl  - RUl)t<TEMPl/FX(l) 

GAMl  - RU2/(XMUXl#TEMP2)  +1.0 
97  IDFF  - 0 

IF  (FFID.EQ.>I)  IDFF  - I 

IF  (FFID.EQ.P7)  IDFF  ■ 2 

IF  (FFID.EQ.P19)  IDFF  • 3 
IF  (FFID.EQ.SP)  IDFF  ■ 4 

IF  (FFID.EQ.CX)  IDFF  - 5 

IF  (IDFF.NE.0)  GO  TO  21 

CALL  PLOT  (27,12,4) 

TYPE  315 

315  FORMAT  (SX,'UNKNQUN  GRAIN  TYPE'/ 

1 5X,* INPUT  TABLE  OF  DISTANCE  BURNT'/ 

2 5X, 'VERSUS  TOTAL  BURN  AREA  OF  CHARGE'/ 

3 5X. 'FIRST  ENTER  NUMBER  OF  POINTS'/ 

4 5X,'TU0  DIGIT  NUMBER  - PLEASE') 

ACCEPT  316,  NPA 

316  FORMAT  (12) 

TYPE  317,  NPA 

317  FORMAT  (5X, ' INPUT', 2X, 12, 2X,' ENTRIES  OF  TABLE'/ 
15X,'DIST  BURNT  - TOTAL  BURN  AREA') 


96 


•I 


0231 

0232 

0233 

0234 

0235 


DO  43  I • l.NPft 

43  ACCEPT  304,  UEBXd),  ABX(I) 

GO  TO  22 

21  TYPE  318,  GL 

310  FORMAT  (5X, 'GRAIN  LENGTH'. F 12. 5) 
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Q236 

ACCEPT  304,  TElf> 

023? 

IF  (TEMP. NE. 0.0)  GL  - TEtl* 

0239 

TYPE  319.  OD 

0240 

319 

FORMAT  (5X.' GRAIN  OUTER  DIAMETER' .F 12. 5) 

0241 

ACCEPT  304.  TEMP 

0242 

IF  (TEMP. HE. 0.8)  OD  ■ TEMP 

0244 

TYPE  320.  PD 

0245 

320 

FORMAT  (5X, 'PERFORATION  DIAMETER' .F12. 5) 

0246 

ACCEPT  304,  TEMP 

024? 

IF  (TEMP. NE. 0.0)  PD  - TEMP 

C 

CALL  PLOT(2?,12.4) 

0249 

TYPE  321.  UID 

0250 

321 

FORMAT  (5X,' INNER  UEB'.F12.5) 

0251 

ACCEPT  304.  TEMP 

0252 

IF  (TEMP. Nl:. 0.0)  UID  • TEMP 

0254 

TYPE  322,  UOD 

0255 

322 

FORMAT  (5X, 'OUTER  UF,B',F12.5) 

0256 

ACCEPT  304.  TEMP 

025? 

IF  (TEMP. NE. 0.0)  UOD  - TEMP 

0259 

22 

TYPE  427.  FRAC 

0260 

427 

FORMAT  (5X,' IGNITION  BRNT  FRACT' .F12.5) 

0261 

ACCEPT  304.  TEMP 

0262 

IF  (TEMP. NE. 0.0)  FRAC  - TEMP 

0264 

CALL  PL0T(2?,12.4) 

0265 

TYPE  887.  PLQ 

0266 

88? 

FORMAT  (5X, 'LOWER  FIT  LIMIT* .F12,5) 

026? 

ACCEPT  304,  TEMP 

0268 

IF  (TEMP. NE. 0.8)  PLO  - TEMP 

02?0 

TYPE  888,  PHI 

02?  1 

888 

FORMAT  (5X,'UPPER  FIT  LIMIT' .F12.S) 

02?2 

ACCEPT  304,  TEMP 

02?3 

IF  (TEMP. NE. 0.0)  PHI  - TEMP 

02?5 

TYPE  531,  UDEV 

02?6 

531 

FORMhT  (5X,'UEB  deviation  iS',Fl2.5) 

02?? 

ACCEPT  304,  TEMP 

02?0 

IF  (TEMP. NE. 0.0)  UDEV  » TEMP 

0288 

TYPE  532,  TDEV 

0281 

532 

FORMAT  (5X,' IGNITION  DEVIATION  IS'.F12.5) 

0282 

ACCEPT  304.  TEMP 

0283 

IF  (TEMP. NE. 0.0)  TDEV  - TEMP 

0285 

SPACEl(10)  - 0.0 

0286 

TYPE  2468.  SPACE  1(10) 

028? 

2468 

FORMAT  (5X.'VENT  AREA  FOR  VENTED  OPERATION' .F12. 5) 

0268 

ACCEPT  2469.  TEW 

0289 

2469 

FORMAT  (E12.0) 

0290 

IF  (TEMP. NE, 0.0)  SPACEKIO)  - TEMP 

0292 

TYPE  5832.  PBLOU 

0293 

5832  FORMAT  (SX,*BLOU  OUT  PRESSURE  LEVEL' .F12.5) 

0294 

ACCEPT  2469.  TEMP 

0295 

IF  (TEMP. NE. 0.0)  PBLOU  - TEMP 

029? 

TYPE  7354,  IBUG 

0298 

7354  FORMAT  (5X,'DIFEQ  TRACE  CONTROL'/ 

98 


0299 

0300 

0301 

0302 


I 5X,*TYPE  FOR  TRACE  0N'»13) 

ACCEPT  ?355.  IBUG 
7355  FORMAT  (ID 
GO  TO  92 
99  CONTINUE 


99 
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0303 

iDPF  • e 

0304 

IF  (FFID.EQ.Pl)  IDFF  • 1 

0366 

IF  (FFID.EQ.P7)  IDFF  • 2 

0300 

IF  <FFID.EQ.P19)  IDFF  - 3 

0310 

IF  (FFID.EQ.SP)  IDFF  • 4 

0312 

IF  (FFID.EQXY)  IDFF  • 5 

0314 

IF  (IDFF.NE.0)  GO  TO  87 

0316 

TYPE  86 

0317 

86 

FORMAT  (5X,'UNi<NQyN  GRAIN  TYPE  IN  STANDARD V 

1 5X,* PROGRAM  HAS  ABORTED') 

0310 

STOP 

0319 

8? 

CONTINUE 

0320 

NPTH  ■ 4 

0321 

TEMPI  - XMUXl 

0322 

TEMP2  - GftMl 

0323 

IF  (GftMl. GE. 1.0)  TEMP2  - 1544. *1 .0/779. /XMJXl/(GflMl  - 

0325 

IF  (XMUXl. LT. 1000.)  TEMPI  ■ F1*XMUX1/(1544.*1.B) 

0327 

P2  ■ 1.0 

0329 

DO  4?  I - 1.4 

0329 

PX(I)  - P2 

0330  .. 

FX(n  ■ Ft  . - . . 

0331 

ETX(I)  ■ ETl 

0332 

XMUX(I)  - TEMPI 

0333 

GAMX(I)  • TEMP2 

0334 

P2  • P2  + 1.0 

0335 

47 

CONTINUE 

3336 

92 

ISKl  - MAXOdSKl.l) 

0337 

IF  (IDFF. EQ, 2. OR. IDFF. EQ. 3)  GO  TO  6138 

0339 

UltlN  • AMINl  (UID.UOD) 

0340 

UQD  - UMIN 

0341 

UID  - UMIN 

11342 

6130 

CONTINUE 

ISK  - 10 

0344 

ISK2  - 5 

0345 

RHO  - RNOC't<.036UU 

0346 

RETURN 

0347 

END 
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SUBROUTINE  HP*.OSS 

C SUBROUTINE  TO  CALCULATE  HEAT  LOSS  NUMBER  - 2 OPTIONS 
COMMON  RUNID(5)>RTITLE(IB),DATEC3).0PERN(5),PR0PN(18), 

1 PSORCEC 18) >PRLOT( 18) ,PREM( 18) ,TIGNR ( 18) ,GAGE( 18) , 

2 SPACE  I (36) , CHGUT. UID, PTEMP, UGHTI , PCORR, BVOL, BTEMP. 

3 SCAL,CALCDN>STIME,FFID,GL,OD.PD,UOD.Fl,DUM,XMUXl,ETl,GAril. 

4 SPACE2 ( 4) , ISK 1 > KTOT> MZ, MY, PMHX, TMAX, NPMA. PPM. DPMAX. IHL. H. 

5 RHOC.RHO,T10.T90,T1690.P9C5).DP9(5).NPTH,PX(20). 

6 FX(20).ETX<20).XMLJX(20).GAMX(20).IDFF.NPA.UEBX(20).ABX(20), 

7 ISK, 1SK2.MEM2.XMI.RP.TZERO.VOL.CP,CST,TM.PIGN,TIGN.PTHEO. 

8 C0Nl,C0N9,MEM3.MEri7.SFAC.L9.P.DP.TST0P, I9. 

9 LS.CONU.TBU.FRAC.Z.SRAT.XTBU 

EQUIVALENCE  (FA ID. SPACE  1 ( 12) ) . (’ ID. SPACE 1(13)). 

1 (XTW.SPACEKU)).  (CAID.SPACEKIS)) 

MEM2  • 1 

CALL  FINDl  (Pr^X.F.PX.FX.NP7H.MF.M2) 

CALL  FINDl  (PMAX.ETA.PX.ETX.NPTH.MEM2) 

CALL  FINDl  <PMAX.TZP.PX,XMLJX.NPTH.MEM2) 

CALL  FINDl  (PMAX.CVP.PX.GAMX.NPTH.MEM2) 

RU  - 1544.*12.n<1.8 
XMU  - RU*TZP/(12.n<F) 

T2ER0  - T2P 
TA  ■ 298 « 

TI  • IZ.H-PCORRi^I/RU 
VOL  • BVQL/(2.54iW<3) 

CA  - 14.7*V0L*29./(298.*fiU) 

Cl  « UGHTI/454. 

CP  • CHGUT/4S4. 

CST  - CA  + Cl 

ETX(20)  - CA/29.  + C^'WI  + CAlD>*fRAC/XMA 
TM  ■ CA>t<290.  CIKTI  + CA ID*FRAC»TAID 
ETX(19)  - CVP*<CST  + FRACXCAID) 

TM  - TM/(CST  FRACHCAID) 

CTOT  - CA  + Cl  CP  + CAID 
RA  - RU/^9. 

RI  - RU/XMl 
RP  - RU/XMLI 
RAID  - RU/XMA 

PTHEO  - <CA*«A*298.  + ClX(RI<«TI  CP*8P*TZER0  CAlDt«AID*TAID) 
I /(VOL  - CTOT*ETA)/1O00. 

TM  - ntKPMAX/PTHEO 
TL  • TI*PMAX/PTHEO 
TLl  - TAID<t<PrtAX/PTHEO 

PIGN  - (CA*TA/29.  + C!*TL/XMI  + FHAC*CAID>t‘TLl/Xr«)»RU/(VOL  - 
I (CST  + FRAC«CAID)>^ETA  - CP/RHO  > ( i . - FRRC)*CRID/RHO)/1000. 
CALL  FINDTP  (PIGN.  TIGN.3.2. 1.2.NPMA.t«M2.t1Z) 

C0N9  - VOL  - CP/RHO  - (1.  - FRAC)>tCAlD/RHO 
IF  (H,EQ.O.0)  GO  TO  10 
IF  (FRAC.EQ.1.0)  RETURN 
10  DPDT  - (PTUEO  - PMAX)/(TMAX  - TIGN) 

H - (VOL  - CTOTil€TA)»CVP#DPDT/RU#l000.*CTOT/ 

1 (ETX(20)  + CP/XMJ  + (1.  - FRAC)*CAID/W1A) 


101 


8841  IF  (IHL.NE.2)  GO  TO  1492 

8043  IF  (FRflC.EQ.l.e)  RETURN 

0045  HTT.!  ■ H 

0046  GO  TO  14S3 

0047  1492  SPACE U36)  - 450. 
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0070 
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0075 
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FORTRAN  IV 

TYPE  372,  S«ICE1<36) 

372  FORMAT  (5X,‘«S€RAGE  B01«  UALL  TEMPERATURE*, FI 2. 5) 

ACCEPT  373,  TElf» 

373  FORMAT  (El2,t) 

IF  (TElP.HE-e.B)  SPACE1C3S)  - TEMP 

H0  - H>t<l080.ySPACE2(2)/((TZERO’4»MAX/1>THEO  + TM)/2. 

I _ SPACE  1(36)) 

UB  ■ (d.  - fBAC)=«:AlD  + CP)/(TMAX  - T1GN)»1000. 

ua  - ije4<«0.a 
H ■ HBAB 
TYPE  932.  HB,  H 
932  FORMAT  (2E16.6) 

SPACE1(34)  -8.0 

1493  IF  (FRAC.EQ.1.0)  GO  TO  68 
TST  - TIGN  - 0.50 

TST  - AMAXHT5T,0.0) 

SUM  - 0.0 
BELT  - 0.02S 
DO  69  I • 1,28 

CALL  FINDTP  (TST,DP.3. 1,3,2. NPMA,MEM2,MZ) 

SUM  - SUM  + iP 
TST  ■ TST  BELT 
69  CONTINUE 

DPB  - SUM*S9. 

PIG  - piGN«teea. 

DYH  - FRAC*CAID/T!GN*l00a. 

DYHL  • DYN 

TSYS  • (C0H9  - (FRACCAID  + CST)*£TA)>*PlG>tTX(20)/RU 
1490  DVN  • DYN-xwa.e 
HTL  - HTLl 

IF  (1HL.EQ.2)  GO  TO  1494 

HTL  • H!*<DVN*6PACE2C2)#(TSYS  - 290.)/'ia00. 

1494  CONTINUE 

DVN  • ((CONS  - (FRAC*CAID  + CST)»£TA)*DP8  + ETX(20)*fiU>»HTL/ 

1 ETX(19))/(ETX(26)i«8U'*<(TAID  - TSYS)/(CST  + FRRC*CAID)  + 

2 RU*TSYS>'XMR  - PIG/'RHO  + P1G»ETA) 

IF  (ABS  (1.  - DYN/DYNL)  - 0.001)  ?IO.310.311 
311  DYNL  • OYN 
GO  TO  1490 
310  CONTINUE 

XNIGN  - 0.8  “■ 

TYPE  1510.  XHIGN 

1510  FORMAT  (SX.'PflUER  ON  IGNITER  FL0U'.F12.S) 

ACCEPT  ISU.  TEMP 

1511  FORMAT  (E12.8) 

IF  (TEr^.NE.0.0)  XNIGN  - TEMP 
SPACE  1(33)  - )«iIGH 
P8  • PI6«>k>«16N 
SPACE  1(34)  - BVN/P0 
TYPE  1496.  SPACE I (34) 

1496  FORMAT  (E16.6) 

68  COllTMUE 


103  I 


104 
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8881  SUBROUTINE  HTFIT 

C SUBROUTINE  TO  FIT  DECAY  TO  FIND  HEAT  LOSS  COEFFICIENT 

0002  COMMON  RUHID(5)>RTITLE(10).DATE(3),OPERN(5KPROPN(ia). 

1 PSORCEC 18:. .PRLOTC 18) .PREMC IB) >TIGNR(  18) .GAGE ( 10) . 

2 SPACE  1(36). CHGUT.UID, PTEMP, UGHTl .PCORR. B'^/OL. BTEl^. 

3 SCAL.CALCON.STIME.FFlD.GL.OD.PD.UOD.Fl.DUM.XMJXl.ETl.GAMl. 

4 SP ACE2 ( 4) . I SK 1 . KTQT. MZ,  MY. PMAX. TMAX. NPMA. PPM. DPMAX. I HL. H. 

5 RHOC.RHO,T10.T90.TlQ90.P9(5).r^S(5).NPTH.PX(20). 

6 FX(20) . ETX(20) .XMUX(20) .GAMX(20) , IDFF, NPA.LJEDX(20) . ABX(20} . 

7 ISK. ISK2.MEM2.XMI.RP.TZER0.VQL,CP,CST.TM.P1GN.TIGN.PTHEQ. 

8 CONI. CONS. MEM3.MEM7.SFA2.L9.P. DP. TSTOP, 19. 

9 L8.C0NU.TBU.FRAC.2.SRAT.XTBU 


8803 

COMMON  /GENE/  UIDP(5).U0DP(5).UCR(5).SFAC(5.5).TDEL(5). 

1 UDEV.TDEV.TUST 

0004 

COMMON  /^LAH/  PLO.  PHI.  ACO.  XNC 

0085 

DIMENSION  Al(500.I).ai(5a0.  D.DKI.  1) 

0006 

20 

r«M2  - I 

0007 

MEK3  • 1 

0008 

NE  - NPMA  + 50/ ISK 1 

0009 

INOT  - (KTOT  - HE)/5a0  * 1 

0010 

CTQT  - CST  + CP 

0011 

RU  « 12.  =*>1544. *1.8 

0012 

CALL  FINDl  (PMAX.  F.  PX.  FX.  NPTH.  rEM2) 

0013 

CALL  FIND!  (PMAX.  TZP,  PX,  XMUX.  NPTH.  r^M2) 

8014 

CALL  FINDl  (PMAX.  CVP.  PX.  GAMX.  NPTH.  11Ef12) 

80  IS 

XMU  • RU»TZP/(12.*F) 

0016 

eOT  • (ETX(28)  ♦ CP/XM1J)*RU 

0017 

44 

J • 0 

8318 

DO  62  1 - HE.KTDT. INOT 

0019 

READ  (2*1)  Tl.  P.  DP 

0020 

CALL  FIND!  (P.  ETA.  PX.  ETX.  HPTX.  MEM2) 

0021 

DETA  • (ETX(MEM2*1)  - ETX(MEM2)  )/(PX(MEM2>t)  - 
1 PX(I€M2)) 

0922 

P • P*1B00. 

8023 

DP  - DP*  1000. 

0024 

DETA  - DETA/1080, 

0825 

T • P>»(VOL  ' CTOT*£TA)/0OT 

0026 

J » J + 1 

9027 

AUJ.l)  T**4 

6828 

01(J.l)  • -DP«(VOL  - CTOTweiA  - CTOTit.DETA>*P)*CVP/RU*lO00 
1 vCTOT/(ETX(20)  ♦ CP/XJU) 

0829 

62 

CONTINUE 

0030 

CALL  SIM€Q  (At.ei.Dl.J. 1. 1) 

0031 

SPACE2{3)  • Dl(l.l)/SPAC£2<2) 

0832 

SUM  • 0.0 

0833 

DO  777  17  - l.J 

0034 

777  SUM  - SUM  ^ BKJ.  I)«e*2 

0035 

RMS  - SQRT(SUM/(J  - D) 

883S 

Ue  - (CP/(TMAX  - TlGN))*108a. 

8037 

U0  - liB*»0.0 

6838 

TAV  - (TZER0**4  ♦ TM**4)/2. 

0039 

TGOQ  - (TZERO  ♦ TM)/2.  - SPACE  1(36) 

lOS 


ae4a 

HI  • 

H«UB>kTGQO  > SP8CE2<3)«TAV 

Q841 

SUMl 

- 8.0 

8042 

SUM2 

- 0.0 

0043 

DELT7  - (TMAX  - TIGN)X10a0. 

0044 

T7  - 

TIGN 
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0045  636 

0046 
8347 
0346 
0349 
0356 
0d52 
0tJ53 

6Q54 

6055  100 

0656 

0057 


CALL  FINDTP<T7,DPT>3, l.3,2*3000*MEM2) 
SUMl  - sum  + DPT*DELT7 
DPT  - AMfiXl (DPT, 0,00001) 

SUM2  • SUM2  + DELT7»DPT#*0.8 
T7  - T7  DELT7 
IF  (T7.LE.TMAX)  GO  TO  636 
FAC7  - SUMl*i«0.8/SUM2 
H - HL/U6/TGQ0 
H • Hn<FAC7 

TYPE  100.  DUl.U,  SPACE2(3),  RMS,  H 
FORMAT  (4E15.5) 

RETURN 

END 
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0001  SUBROUTINE  SUMUP 

C SUBROUTINE  TO  PRINT  SUMMARY  SHEET  OF  ANALYSIS 

0002  COMMON  RUNIO(5URTITLEn8)>DATE(3)AOPERN(5),PROPN(l0). 

1 PSORCEdB)  .PRLOT(  18)  ,PREM(  18)  ,TIGNR(  18)  .GAGE(  18) . 

2 SPACE  I <3G) , CHGUT, U ID, PTEMP, UGHT I , PCORR. BVOL , BTEMP, 

3 SCAL,CALC0N,STIME,FFID,GL,0D,PD,U0D,FI,DUM,XMUX1,ET1.GRM1, 

4 SPACE2 (4) , ISKi , KTOT, M2, MY, PMAX,TMAX, NPMA, PPM, DPMAX, IHL,H, 

5 RHQC.RHO,T10,T3Q,T1090,P9(5),DP9(S),NPTH.PX(20), 

6 FXC20)  ,ETX(20) , XMUX(20) ,GAMXC20) , IDFF,NPA, UEBX(20) , A0X(20) , 

7 ISK, ISK2,MEM2,XMI,RP,TZER0,V0L,CP.CST,TM,PIGN,TIGN,PTHE0, 

8 CaNl,C0N9,MEM3,MEM7,SFAC,L9,P,DP,TST0P, 19, 

9 L0,CQNU,TBU,FRAC,2,SRAT,XTBU 

0003  DIMENSION  ABC  (5) 

0004  DATA  ABC/'CONS','TANr.'PROP',*ORTr,'ONALV 

0005  PRINT  100 

3006  103  FORMAT  (IHl///) 

3007  PRINT  101,  CRUNID(I),  I - 1,5) 

0000  101  FORMAT  (10X,'RUN  IDs*23X,5A4) 

0003  PRINT  102,  <RTITLE(I),I  - 1,18) 

0010  102  FORriAT  CiaX,' RUN  TITLE :",2SX,  18A4) 

0011  PRINT  103.  (DATEd),  I - 1,3) 

8012  103  FORMAT  ( ISX, 'DATE:* ,30X,3A4) 

3013  PRINT  134,  (OPERN(I),  1 - 1,5) 

0314  104  format  a 0X.' operator !',26X,SA4,') 

0315  PRINT  105 

0016  105  FORMAT  (I 0X.' PROPELLANT  DATA'* /) 

001?  PRINT  136.  (PROPNd),  I • I.  IB) 

0010  186  FORMAT  (IBX. 'TYPE5'.33X, 10A4) 

8019  PRINT  10?,  CHGUT 

0030  107  FORMAT  d0X,'UEIGHT  (GK6) : ' ,22X.F12.5) 

0021  PRINT  100,  RHOC 

8022  108  FORMAT  d0X. 'DENSITY  (GM/^^C) 19X,F12.S) 

0023  PRINT  109.  PTEMP 

0024  109  FORMAT  d0X.' INITIAL  TEMPERATURE  (DEu  to  !',?X,F12, 5) 

0025  PRINT  110,  'PRLOTd),  I • 1.18) 

0026  110  FORMAT  d0>C'LOTj'.31X,  1BA4) 

002?  PRINT  111.  <PSORCE<I),  I - 1.10) 

0026  III  FORMAT  dOX. 'SOURCE :',28X.  18A4) 

0029  PRINT  112.  FFID 

0030  112  FORMAT  d0X. 'GRAIN  TYPE:',24X.A4) 

0031  PRINT  113.  GL.  OD.  PD 

9032  113  FORMAT  d9X. 'LENGTH. OD, ID  (IN) , 17X,3(F12.S.',')) 

0033  PRINT  114.  UID,  UQD 

0034  114  FORMAT  (10X.' INNER  UE8.0UTER  UE8  dN) 10X.2(F12,5.'.')) 

0035  PRINT  I IS.  FXd) 

0036  115  FORMAT  d0X, ' THEORETICAL  IMPETUS  <FT-LB/L0) s'.4X.FI2,5) 

0037  PRINT  ns.  TZERO 

0030  116  FORMAT  (lOX.'FLAfE  TErt’ERATUREs', 17X.F12.S) 

0B39  PRINT  117.  XMUXl 

0040  117  FORMAT  d0X, 'AVERAGE  MOLECULAR  UElGHT  OF  PROD: ' .2X,F 12.5) 

0041  PRINT  ne.  ETXCl) 

0042  lie  FORMAT  ( 10X. 'CO-VOLUPE  (CU  IN/LB) J*. 14X,F12.5> 
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0043  PRINT  U3,  Cflm 

0044  119  FORMAT  U0X. ‘GAMMA  (RATIO  OF  SP  NTS) UX,F  12.5) 

0845  PRINT  120.  <PREM(I),  I - 1,18) 

0046  120  FORMAT  < IBX.* REMARKS t'.27X.10A4/) 

8847  PRINT  121 
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8048 

0849 

0050 

0051 

0052 

0053 
8054 

0055 

0056 

0057 

0058 

0059 

0060 
0061 
0062 

0063 

0064 

0065 

0066 

0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 

0078 
0379 
0080 
6381 
0082 
0083 
8084 

0085 

0086 

0087 

0088 
0689 
0030 

0091 

0092 

0093 

0094 

0095 
6036 

6097 

6098 


121  FORMAT  <10X,' IGNITER  DATAiV) 

PRINT  122.  (TIGNR(I),  I ■ UI8) 

122  FORMAT  (10X.*TYPE:'.30X. 10A4) 

PRINT  123.  (JGHTI 

123  FORMAT  (10X.'UEIGHT  (GMS) ;*.22X.F12.5) 

PRINT  124.  PCORR 

124  FORMAT  ( 10X. ' IMPETUS  (FT-L0/LB)  J M6X,F12.5/) 

PRINT  125 

125  FORMAT  (10X. 'EQUIPMENT  DATA*/) 

PRINT  126.  BVOL 

126  FORMAT  (IBX.'BOre  VOLUME  (CC) i ' . 18X.F12.5) 

PRINT  127.  BTEMP 

127  FORMAT  <10X.'BDMB  TEMP  (DEG  K) s ' 17X,F12.5) 

PRINT  128.  (GAGEd).  1 - 1.18) 

128  FORMAT  (10X, 'GAUGE  TYPE: ' ,24X.  18A4) 

PRINT  t29/SCAL 

123  FORMAT  ( 10X. 'CALIBRATION  FACTOR  (PC/PSI) :'.7X.F12.5/) 

. PRINT  130 

130  FORMAT  (18X.' RESULTS:'/) 

PRINT  131.  PTHEO 

131  FORMAT  <10X.' THEORETICAL  MAX  PRESS  (KPSIA) :'.5X.F12.5) 

PRINT  132.  PMAX 

132  FORMAT  ( 10X. 'OBSERVED  MAX  PRESS  (KPSIA):'.  0X.F12.5) 

PRINT  133.  PIGN 

133  FORMAT  ( 10X.' IGNITER  PRESSURE  (KPSIA) 10X.F12.5) 

PRINT  134 

134  FORMAT  (/10X.' IGNITION  TIME  INFORMATION:'  ) 

PRINT  135.  T10 

135  FORWT  (10X.'Tir€  TO  10X  PMAX  (MSEC) 11X.F12.S) 

PRINT  136.  T90 

136  FORMAT  (10X.'TIME  TO  90X  PMAX  (MSEC) : ' . 1 1X.F12.5) 

PRINT  137.  TMAX 

137  FORMAT  (l0X.'TiriE  TO  I08x  PMAX  (MSEC) : ' . 10X.F12.5) 

PRINT  130.  T1090 

130  FORhWT  (19K.'TIME  FROM  I8x  TO  90X  PMAX  (MSEC)  :'.2X.F12.5> 
PRINT  140 

140  FORINT  (IHI////10X.  'QUICKNESS  INFORMATION: */) 

PRINT  141.  DP9(1) 

141  FORMAT  (10X.'POOT  AT  .250  PMAX:' .5X.F12.S) 

PRINT  142.  DPS(2) 

142  FORMAT  (10X.'PDOT  AT  .375  PMAX: ' ,5X,F 12.5) 

PRINT  143.  DPgO) 

143  FORrWT  (10X.'PDOT  AT  .500  PMAX:-* .5X,F12.5) 

PRINT  144,  DP9(4) 

144  FORMAT  (10X,’PDOT  AT  .625  PMAX: ' ,5X.F12.5) 

PRINT  145.  CP9(5) 

145  FORMAT  ( 1 0X. ' AVERAGE  PDOT : ' . 1 0X. F 1 2 . 5/) 

PRINT  139.  DPMAX.  PPM 

139  FORMAT  ( 10X.' MAXIMUM  PDOT  (rt»SI/SEC> :' .FU .5.5X. 'OBSERVED  AT'. 
1 ' P • (KPSIA):'.F12.S/) 

GO  TO  (1.2).  IHL 
2 PRINT  146.  (ABCd).  I • 1.2) 
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0893  146  FORMAT  C “ HEAT  LOSS  OPTION i S 2X*  ZA4) 

eies  GO  TO  3 

eiei  l PRINT  147.  (ABC  (I).  I - 3.5) 

0102  147  FORMAT  (10X.*HEAT  LOSS  OPTIONi%2X.3A4) 

0103  3 PRINT  148.  H 
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eeei  subroutine  reduce 

C DRI'^R  FOR  DIFFERENTIAL  EQUATION  SOLVER 

0082  COmON  RUHID(S>.RTITLE(10)«DATE(3).OPERN(5).PROPN(18). 

1 PSORCE( 18) .PRLOTC 10) ,PREM( !8) «TIGNR( 16) ,GAGE ( IB) . 

2 SPACE! (36). CHGUT,UID,PTE^P.UGHTI,PCORR.BVOL.BTEMP. 

3 SCAL. CALCON. STIME. FF ID. GL. OD>  PD. UOD. F 1 . DUli. XMUXl . ETl . GAHl . 

4 SPACE2  ( 4) . I3K I . KTOT. MZ. W.  PMAX.  TMAX.  NPflfl.  PPM.  DPMAX.  IHL. H. 

5 RHOC.RHO.T10.T90.T1090.P9(S).DPS(5).NPTH.PX(20)- 

6 FX(20) .ETX(20) .XMyX(2e) .GArt<(20). IDFF.NPA.UEBX(20) . ABX(20) , 

7 ISK. ISK2.rCti2.XMI.RP.TZER0.V0L.CP.CST.TM.PlGN.TIGN.PTHE0. 

8 CONI. C0N9. hEM3. MEt17. SFAZ. L9 . P. DP. TSTOP.  19. 

9 L8.C0HU.T6U.FRAC.Z.SRAT.XTBU 

0803  COftIDN  /GENE/  UIDP(S).iJ0DP(S).UCR(S).SFAC(5.5).TDEL(5), 

I UDEV.TDEV.TUST 
0084  DIftNSION  TPRNT(2) 

0005  TPRNT(l)  — STIME«l.0E-03 

6006  TPRNT(?)  - 10000. 

000?  TSTOP  - TMAX 

0008  IF  (SPACEI(10).NE.0.0)  TMAX  • 290a>*STIf€*l .BE-03 

0010  IF  (SPACEU10).NE.0.8)  H • 0.0 

0012  JDO  - 5 

0013  IF  (TDEV.EQ.8.0)  JDO  • 1 

0015  JDO  • JDO  + 7 

0016  CALL  EVES  (JDO.  TPRNT) 

8017  ENDFILE  2 

0018  ENDFILE  6 

0019  RETURN 

0020  END 


113 


RT-U 

0081 

0002 


0003 

0004 

0005 

0006 

0007 

6008 

0009 

0010 
0011 
0012 

0013 

0014 

0015 

0016 

0017 

0018 
0019 
0021 
0023 
8024 

j827 

0028 

0829 

0030 

0031 

0032 

0033 

0034 

0035 
@037 

0038 

0039 
8040 
0842 
0043 
0644 

0045 

0046 
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SUBROUTINE  SETUP  (T.  Y»  BIG.  N) 

COMMON  RUNID(5) .RTITLE( 10) ,DATE(3) ,0PERN<5) .PROPNI 10) , 

1 PSORCEI IB) ,PRLOT( 18) .PREMC IB) .TIGNRI 18) ,GRGE( 18) , 

2 SPACE  1(36). CHGUT, UI D. PTEMP. UGHT I . PCORR. 8V0L . BTEMP, 

3 SCAL,CALC0N.STIME.FFID.GL.0D.PD.LI0D,F1.DUM.XML)XI.ETI,GAM1, 

4 SPACE2  ( 4) . ISK 1 . KTOT.  MZ.  MY.  PMAX.  TflAX.  NPMA , PPM. DPMAX.  IHL.  H, 

5 RHOC.RHO.T10.T90.T1Q90.P9(S).CP9(5),NPTH.PX(20).. 

6 FX(20) .ETX(20) .XMIJX(2B) .GAMX(20) , IDFF.NPA.UEBX(2B) . ABX(20) . 

7 ISK. ISK2.MEM2.XMI.RP.TZERQ.V0L,CP.CST.TM.PIGN.TIGN.PTHE0. 

8 CQN1.C0N9.MEM3.M£M7.SFAZ,L9.P,DP,TST0P.  19. 

9 L9 . CONU. TBU. FR AC . Z , SRAT, XTBU. RATE . RL 

COrtlQN  /GENE/  UIDP(5).UODP(5).UCR(5).SFAC(S.5),TDEL(5). 

1 IJDEV.TDEV.TUST 
COMMON  /KEVIN/  TE(35).  XT(35) 

DIMENSION  JD(2) 

EQU IVALENCE  ( JD  ( I ) .SPACEH35) ) . (CA ID. SPACEl  ( 15) ) . 

1 (CPGM2,  SPACE  U 16) ).  (P8L0U.  SPACE  KID) 

DIMENSION  T(2).  Y(2).  SIG(12.3) 

Dlf^NSlON  CPR0(5).PRQ(5).DDU(5) 

DATA  PRO/.  1.  .2.  .4.  .2. . lADDUz-1 . 17741. -.83255.  .0.  .83255. 1 . 17741/ 

MEM2  - I 

l«rt3  » I 

Cl  • IJGHTI/4S4. 

CALL  FlNDl  (PMAX.F,PX.FX,NPTH.MEri2) 

CPGM2  - GAMl*(2./(GAMl  + I . ) )>Mt(  (GAMl  -f  l.)/(GAMl  - I.)) 

A8X(20)  - l.aE+05 
MEM7  • 1 
JD(I)  • 5 
JD(2)  - S 

IF  (TDEV.EQ.0.0)  JDd)  • 1 

IF  (UDEV,EQ.0.0)  J0(2)  - 1 

C0H9  - VOL  - (CH6UT/454.  ♦ CAID)/RH0 

IF  (ID^F.EQ.0)  GO  TO  10 

19-0 

S2ERQ  - 0.0 

SUM9  - 0.0 

JDO  • JD(l) 

DO  75  J • l.JDO 
Jl  - J + 7 
SIG  (J1.2)  - 0.001 
SIG  (J1.3)  • l.0E->e5 
CP ART  - CP»PRO(J) 

IF  (JDO.EO.l)  CPART  - CP 
IDO  • JDC2) 

DO  7S  I - 1.  IDO 
CPRO(I)  - CPART>i4’R0(l) 

IF  (IDO.EQ.D  CPRO(I)  • CPART 
LilDP(I)  - UID  + DDUtDHJDEV 
UODP(l)  - UOD  - DDU(l)+yOEV 
19-0 

GO  TO  (1.2. 3. 4.5).  IDFF 

1 CALI  FORM1(0.0.  I9.OD.PD.U1DP(I).UODP(I).V0,GL.S.UCR(I))  ' 
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0047 

0040 

0049 

0050 

0051 


GO  TO  39 

2 C0LL  PORM7(0.8.I9,OD<PD.UI9P(n.UQSP(l}.VO.GL<8.UCR(I)) 
GO  TO  39 

3 CALL  P0RMl9(e.e>I9«01)>PD,UI9P(I).U0DP(I)<Ve.GL.8«UCR(I)) 
GO  TO  39 
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0052  4 CALL  FORrtSP(0.0.I9.OD,PD>UIDP(I).UODP(I).V0.GL.S.UCR(I)) 

0053  60  Tt)  39 

0054  5 CALL  FORhCV(0.0* I9-0D>PD.UIDP(I) ,UODP( 1) .VB,GL,S,UCR(I)) 

0055  39  SFAC(J,I)  - CPRO(n^(RHO*V0) 

0056  SUMS  - SUMS  SFAC(J.I) 

0057  S2ERQ  - S2ER0  + SFAC(J.I)*S 

0050  SI6(J1.3)  - AMINt  (SIG(JU3).  UCR(I)) 

0059  76  CONTINUE 

0060  TDEL<J)  ■ DDU(J)#TDEV 

0061  IF  (J.EQ.n  TDS  - TDEL(J) 

0063  TOEL<J)  • TDEL(J)  - TDS 

0064  75  CONTINUE 

0065  10  CONTINUE 

0066  IF  (IDU.EQ.l)  GO  TO  723 

0060  DO  723  I - U4 

0069  K • I + I 

0070  DO  723  J • K.S 

0071  IF  (UCk(J).GE.UCR(I))  GO  TO  723 

0073  A1  •>  UCR(I) 

0074  A2  • UlDPd) 

0075  A3  - UODP(I) 

0076  DO  724  L - US 

0077  A4-SFAC(L.n 

0070  SFACCL.l)  • SFAC<L.J) 

0079  SFAC  (L,J)  • A4 

6080  724  CONTINUE 

0081  UCRCI?  • UCR(J) 

0092  UlDPd)  • UIDP(J) 

0093  L*0DP<1)  • UODPCJ) 

0084  UCj?<J)  - Al 

08B5  UIDP(J)  • A2 

0086  l!ODP(J)  • A3 

0907  723  CONTINUE 

0038  L3  - 0 

0089  Td)  - TIGN 

0030  IF  dDFF.EQ.9)  UCRd)  - UEBX(NPA) 

0092  IF  (IDFF.EQ.a)  UCR(2)  - UE8X(NPA) 

0894  IF  dDFF.EQ.8)  UCR<3)  • IJE8X(NPA) 

0098  IF  dDFF.EQ.0)  UCR(4)  - LE8X(NPft) 

0098  IF  dDFF.EQ.0)  UCRtS)  - IJEBXINPA) 

0109  IF  dDFF.EQ.0)  S2ER0  «*  AaXCl) 

0102  NUm  - SUt19 

0103  IF  dDFF.EQ.0)  NUfS  - \ 

0103  PRINT  160.  S2ER0.  NUMB 

0105  160  FORmT  (/lOX.' INITIAL  SURFACE  AREA  (SO  IN);*.FI2.5/ 

1 iOX. 'NUMBER  OF  GRAIKS.’Mex.lB) 

0107  Yd)  • CST  - UGHTI/^. 

0106  Y<2)  - yGHTl/454. 

0109  Y(3)  • FRAC*CAID 

0110  Y(4)  • 0.6 

01 U Y<5)  • Y(3) 

01 17  Y(6)  • 0.0 
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aii3 

eu4 

6115 

0116 

0117 


Y(7)  • PIGN 
SIG(U2)  ■ 
SlGfS.?)  • 0.080i 
S1G<3»2)  • 0.0001 
SJ5<4,2)  • 0.081 
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6118 

6119 

0128 

8121 

0122 

0123 

0125 

0126 
ni27 
0128 

0129 

0130 

0131 

0132 

0133 
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SIG(5,2}  - 0.0001 
SIG(5.3)  - CAID 
SIG(6>2)  - 6.001 
SIG(7.3)  - PBLOU 
SPACE1(31)  - 0.0 
IF  (1HL.EQ.2)  RETURN 

DN  - (l>ShTI/'454.  + FRAC*CAID)/TIGN*1080. 
G ■ H*DN**0.0 

CALL  UALTEM  (G.T(l)>299..'m) 

CALL  UALTEM  (G.  T( 1).290.,TM) 

SFA2  - TM 
DO  1729  17  - U35 
1729  TE(17)  - XT(1?) 

RETURN 

END 
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8001 

8002 


0003 

0004 
0085 

0006 

0007 

6008 

0609 

0010 

0011 

0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 
0021 
0022 

0023 

0024 

0025 

0026 

0027 

0028 
0029 
0830 

0031 

0032 

0033 
0834 
0035 
0'036 
003? 
0039 
9040 
0941 
0042 


SUBROUTINE  DIFEQ  (T,  Y,  DY,  N,  TPR) 

COmON  RUNIDC5),RTITLE(18).DflTE(3),0PERN(5).PR0PN(18), 

1 PSORCEC 18) ,PRLOT( 18) ,PREM( 18) ,TIGNR( 10) ,GflGE< 18) , 

2 SPACE  1(36), CHGUT, DID. PTElf* , UGHT I , PCORR , SVOL , BTEMP , 

3 SCAL,CRLC0N,STINE,FFID,GL.0D,PD,U0D,Fl,DUH,XMUXl,ETl,GflI1U 

4 SPftCE2 ( 4) > ISK 1 , KTOT, M2, MY, PMAX, TMAX, NPMA , PPM, DPMAX, IHL , H, 

5 RHOC,RHQ,T10,T90,Tia90,P9(5),DP9(S),NPTH,PX(29), 

6 FX(20)  ,ETX(20)  ,XML1X(20)  ,GflMXC20) , IDFF,NPfi,b£BX(20)  ,ftBX(20) , 

7 ISK, ISK2,MEM2,XMI,RP,TZER0,V0L,CP,CST,TM,PIGN,TIGN,PTHE0, 

0 C0N1,C0N9,MEM3,MEM7,SFA2,L9,P,DP,TST0P,  19, 

9 L8 , CONU, TSU, FRAC , 2 , SRAT, XTBU, RATE , RL 

EQUIVALENCE  < IMUD( 1) ,SPACE2(4) ) , ( IBUG, IMUD(2) ) 

DIMENSION  IMUD(2) 

COMMON  /GENE/  UIDP(5),U0DP(5),UCR(5),SFAC(5,5),TDEL(5), 

1 UDEV,TDEV,TUST 

COMMON  /KEVIN/  TE(35),  XT(35) 

DIMENSION  T(2),  Y(2),  DY(2),TPR(2) 

•DIMENSION  ARATJ(S) 

DIMENSION  JD(2) 

EQU IVALENCE  ( JD ( 1) ,SPACE 1 (35) ) , (TA ID, SPACE  1(13)), 

1 (XMA,SPACE1 ( 14) ) , (CPGM2,SPACE 1(16)) 

CALL  F INDTP  (T( 1) .P,3, 1,2,2, 3000.MEM3,  M2) 

CALL  FINDTP  (T(1),DP,3, 1,3.2,3000,MEM3,  M2) 

CALL  FINDl  (P,  ETA,  PX,  ETX,  NPTH,  MEM2) 

CALL  FINDl  (P,  F.  PX,  FX,  NPTH,  MEM2) 

CALL  FINDl  (P,  T2P,  PX,  XMUX,  NPTH,  MEM2) 

CALL  FINDl  (P,  CVP,  PX,  GAMX,  NPTH,  MEK2) 

RU  - 1544.>t<1.0 
TI  - PCORR»XMI/RU 
XIU  - RU#T2P/F 
JDO  • JD(1) 

IDO  - JD(2) 

DETA  - (ETX(MEM2+I)  - ETX(MEM2) )/(PX(t1EM2+l)  - PX(."EM2)) 

DIMP  - (FX(MEM2+1)  - FX(MEM2) ) /(PX(MEM2+1)  - PX(M£M2)) 

DCVP  - (GAMX(MEM2+l)  - GhMX(MF.M2)  )/(PX(MEM2+l)  - PX(MEM2)) 

DTZP  - (XMUX(MEM2U)  - XMUX(MEM2)  )/(PX(MEM2+l)  - PX(r€M2)) 

DXM*  - RU)t<(DT2P/F  - T2P>*<DIMP/(F>M(2)) 

P - P*1000. 

OP  • DPHiaOB. 

XTl  - 10000. 

DYl  - DY(S)  + DY(6) 

DN  • AMAXl  (DYI>HI000..  0.000601) 

DETA  - DETA/iaoe.x^DP 

DIMP  • DIMP/1000.*DP 

DCVP  - DCVP/1000.>t‘DP 

DT2P  - DT2P/1600.'«DP 

DXMP  - DXMP/i000.*DP 

IF  (IDFF.NE.0)  GO  TO  10 

CALL  FINDl  (Y(8).AB,UEBX,ABX,NPA,PEM7) 

30  TO  40 
10  AAB  - 6.0 
UIP  • 8.0 
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0043  DELJ  - TIGN 

0044  TDLJ  -0.0 

0045  . DELI  - ABX(2e} 

0046  DO  75  J - l.JDO 

0047  Jl  • J + 7 
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9049 

0049 

0050 
0052 
0034 
0055 
3056 
0057 

0059 
G0S9 

0060 
0061 
0062 
0363 
0064 
0066 

0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 
0070 

0079 

0080 
0001 
0082 

0003 

0004 

0085 

0086 
0087 
0080 
0009 

0090 

0091 

0092 

0093 
0894 

0095 

0096 
0090 

0099 

0100 
0101 
0182 
0103 


ARATJ<J)  - 1.0 

ARATS  » 1.0 

IF  (J.EQ.l)  GO  TO  1970 

IF  (TnEKJl.EQ.TDLJ)  GO  TO  l970 

ARATJ<J)  - (T(l)  - DELJ)/(TDEL(J)  - TDLJ) 

ARATJ(J)  - AMINl  (AliAls:i(ARATJ(J),0.0),  1.0) 

ARATI  ' (T(t)  - DELI)^(TDEL(J)  - TDLJ) 

ARATI  - AMINl  (ftMAXUARATI,0.0),  1 .0) 

ARATI  - 1.0  - ARATI 
1978  CONTINUE 
AO  - 0.8 
ARATX  - 1.0 
UCRI  - UCR(I) 

DO  76  I ■ l.IDO 

IF  (UCR(l).EQ.UCRI)  GO  TO  486 

ARATK  - (Y(Jl)  - UCRD/CUCRd)  - UCRI) 

ARATK  - AMINl  (AMAXl  (ARATK.  0.0).  1.0) 

ARATK  - 1,0  - ARATK 
486  CONTINUE 

GO  TO  (I, 2.3. 4.5).  IDFF 

1 CALL  F0RMl(Y(Jl).I9.0D.PD.iJIDP(n.U0DP(I).VB.GL.S,UCI<(I)) 

GO  TO  39 

2 CALL  FORM7(Y(Jl).I9.OD.PD,UIDP(I).UODP(I).V0.GL,S,UCR<I)) 

GO  TO  39 

3 CALL  F0RM19  (Y(  Jl) . 19.0D.PD,UIDP(  I)  .l.X)DP(  I) , V0.GL.S.UCR(  D) 
GO  TO  39 

4 CALL  FORMSP  (Y( Jl) . I9.QD.PD,UIDP( I) ,UODP( I ) .V0.GL.3.UCR( D) 

GO  TO  39 

5 CALL  FORfCY  (Y(Jl). 13,OD.PD.UlDP(I).yODP(l).V0.GL.S,UCR(I)) 

39  CONTINUE 

AB  - AB  + SFAC(J.I)'i<S*»<ARATK 
UCRI  - UCR(l) 

76  CONTINUE 

AAB  - AAB  + A8'«4lRATJ(J)'fARATi 
DELJ  - TIGN  + TDEL(J) 

DELI  • ABX(20)  * TDELIJ) 

TDLJ  - TDELCJ) 

75  CONTINUE 

40  CONTINUE 

US'<‘S  « Y(t)  + Y(2)  ♦ Y(3)  + Y(4) 

VSYS  - CQN9  + (Y(5)  + Y(6))/RH0  - USYS*€TA 

TQTMOL  - Y(l)/29.  + Y(2)/XMI  + y(3)/XI«  + Y(4)Am 

RU  - 1544. *1.8 

RSYS  - PU*T0T110L/USYS 

TSYS  • P*VSYS/<12.»RSYS»<USYS) 

IF  (IBUG.EQ.B)  GO  TO  7354 
TYPE  2936.  TSYS. P. VSYS, RSYS.LiSYS 
2936  FORMAT  (SE14,5) 

7354  CONTINUE 

CSTAR  • SQRT(32.174M?SYS*TSYS/CPGM2) 

DHN  ■ 32.  l74*P>*SPACEU3l)/CSTAR/'i00a. 

GO  TO  (20.21).  IKL 


0104  20  DO  1739  17  - U3S 

0105  ir39  XT(17)  - TE<I7) 

0106  G • H*DHno*<0.a 

0197  CALL  UflLTEM  (G.T<2),SFAZ,TSV8> 

9108  97  G ■ H*DN**0.8 
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0109 

0110 

0111 

0112 

0113 

0114 

0115 

0116 
0117 

one 

0119 

0120 


0121 

6122 

0123 

0124 

0125 
0127 

0129 

0130 

0131 

0133 

0134 

0135 
0137 

0139 

0140 

0142 

0143 

0145 

0146 

0147 
0(48 

0149 

0150 

0151 

0152 

0153 

0154 

0155 

0156 

0157 
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ARP  - G>t<<TSYS  - Kr(l))/'1000.>xSPACE2(2) 

ARP  - ARP  + SPACE2<3)!*‘(TSyS»*4  - XT( l)iM<4)/1000.# 

1  SPACE2(2) 

GO  TO  22 

21  ARP  • H 

22  CONTINUE 

UYd)  - -Y(l)/USYS>t<DNN 
DY(2)  - -Y(2)/USYSi*<DNN 
DY(5)  - SPACE1(34)*P>m^PACE1C33) 

DY(3)  - DY(5)  - Y(3)AISYS*DNN 
RUT  - DY(1)  + DY(2)  + DY(3)  - Y(4)/USYS*DNN 
RMOL  - DY(l)/29.  + DY(2)/^I  + DY(3)/XMA  - Y(4)/USYS*DNN/mJ 
I -Yi:4)n<DXMP/'(XMU»>K2) 

DY(6)  • CTSYS>kRUT  - TAIDH^DYCS)  + USYSn<TSYS*DCVP/CVP 

1 + ARP/CVP  + GAM1*TSYS*DNN  + P/(12.*RSYS) 

2 *<DY(:5)/RH0  - USYS+DETA  - ET0>4RLJT  + VSYS*DP/P  - VSYS/ 

3 USYS-«(RU/RSYS*(RMOL 

4 - TQTH0L'!<RUT/tdSYS))))/'(T2P  - TSYS  + P/(12.*RSYS)w 

5 (ETA  - d.-RHO  + VSYSAJSYS*(  1 . + RU/1?SYS*(  1 ./mj  - TOTMQL/ 

6 USYS)))) 

DYl  ■ DY(5)  + DY(6) 

DY(4)  ■ DY(6)  - Y(4)/IJSYS*DNN 

DY(7)  - DP 

GO  TO  (92.93).  IHL 

92  IF  (A8S(1.  - DYl>t<l000./DN).LE. 0.001)  GO  TO  94 
IF  (DYl. LE. 0.8)  GO  TO  93 
DN  - DYl*1000. 

GO  TO  97 

94  IF  (ABSd.  ~ XTd)/)Cn).LE. 0.0001)  GO  TO  93 
XTl  - XTd) 

GO  TO  20 

93  IF  (Y(8).EQ.0.0)  S7ER0  - AA0 
IF  (Y(8) .EQ.0.0)  ZZERO  - Y(6) 

SRAT  - AA0/SZERO 
IF  (TDEV.NE.0.8)  SPAT  - SRAT^'IO. 

2 ^(Y(6)  - ZZERO) *454. /TIHGUT 
IF  (AA8.GT.0.0)  GO  TO  675 
RATE  - RL 
GO  TO  676 

675  RATE  • DY(6)/(RH0*«A8) 

676  CONTINUE 

DO  745  I - l.JDO 
J • I * 7 

DY(J)  - RATE*ARATJd) 

745  CONTINUE 

Di  ■ 

RATE  - RATE*1030. 

FX(20)  - TSYS 

RETURN 

END 
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6002 
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0007 
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0009 

0010 
0012 
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0017 

0018 
0020 
0021 
0022 

0023 

0024 
0026 

0027 

0028 
0030 
0032 
0013 

0035 

0036 

0037 

0039 

0040 

0042 

0043 

0044 

0045 

0046 

0047 
0049 
0051 
0952 
0053 
0055 


SUBROUTINE  PRINT  <T,  Y,  DY,  N,  TPR) 

COMMON  RUN  ID (5) . RTITLE < IB) . DATE (3) . OPERN (5) , PROPN ( 19) . 

1 PSQRCEUa)  ,PRLOT(  IB)  ,PREM(  18)  .TIGNR(  18)  ,GAGE(  10)* 

2 SPACE1(36),CHGUT.UID*PTEMP.UGHTI>PC0RR.BV0L*BTEMP* 

3 SCAL.CALCON*STIME.FFID*GL>OD.PD.UOD*FI>DUM.XMJX1.ET1*GAM1* 

4 SPACE2 (4) , ISK 1 , KTOT, MZ,  MY, PMAX, TMAX, NPMA , PPM, DPMAX, IHL, H, 

5 RHQC,RNQ,T10,T90,T109O,P9(S),DP9(5),NPTH,PX(20), 

6 FX(20> ,ETX(20) ,XMUX(20) ,GAMX(20) , IDFF,NPA,LJEBX(20) , A0X(20) , 

7 ISK.  ISK2,MEM2,XMI,RP.T3ERO,VOL,CP,CST,TM,PIGN,TIGN,PTHEO, 

8 CONI. C0N9 , MEM3. MEM7. SFA2, L9. P, DP, TSTOP, 19, 

9 L8 > CQNU, TBU. FRAC, Z, SRAT. XTBU,  RATE , RL 

COMMON  /GENE/  UIDP(5),U0DP(5),UCR(S),SFAC(5,5),TDEL<5), 

I UDEV.TDEV.TUST 
COmON  /KEVIN/  TE(35),  XT(35) 

DIMENSION  T(2),  Y(2),  DY(2),  TPR(2) 

DIMENSION  JD(2) 

EQUIVALENCE  (JD(l).SPACEl(35)) 

DATA  XYZ/1010Q0./ 

JDO  - JD(l) 

IF  (L9,EQ.0)  L6  - 0 

IF  (L9.EQ.0)  L8  • 0 

IF  (L9.EQ.0)  L7  - 0 

DO  1729  17  - 1.3S 

1729  TE<I7)  • XT(I7) 

IF  (TPR(3)  .EQ.5.0)  SPACEl(34)  • 0.0 
SFAZ  • FX(20) 

I5TP  • 0 

DO  669  I - I, JDO 
J - I + 7 

IF  (y(J).LT.UCR(I))  GO  TO  669 
ISTP  - ISTP  I 
669  CONTINUE 

IF  ( ISTP. EQ. JDO)  N - 8 
IF  (L9.EQ,0)  GO  TO  777 
>3)0  « JDO  ♦ 7 

IF  (TPR(3) .LT.8.0.OR.TPR(3).GT.XDO)  GO  TO  4336 
ABX(2Q)  • TM) 

4936  CONTINUE 

XDO  - XDO  + 1.0 

IF  (TPR<3)  .Ea.7.0)  SPACEKSl)  - SPACE  I (10) 

IF  (TPR(3) .HE.XDO)  RETURN 
777  CONTINUE 
P - p/iooa. 

DP  • DP/1080. 

DIST  • 2.«»cY(8) 

DPRAT  - DP/PtlAX 

IF  (M0D(L9. ISK).NE.0)  GO  TO  676 
IF  (riOD(L7.40).NE.0)  GO  TO  677 
PRINT  100.  (RUNlD(l),  1 - l.S) 

677  17  - L7  + I 

IF  (riOD(L6.30).NE.0)  GO  TO  832 
CALL  PLOT  (27. 12,4) 
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6056 

832 

U6  - L6  + 1 

0057 

PRINT  10U  T(l).  P,  DP>  RATE.  Z.  SRAT.  DIST,  DPRflT 

0050 

TtFE  302.  T(l).  SFA2.  XT(1) 

0059 

302 

FORMAT  (5X,  3E15.6) 

0060 

676 

CONTINUE 
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0061 

0063 

0064 
y'  0065 

0066 

0067 

0068 
0069 
0071 

0073 

0074 


0075 

0076 
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IF '(tioDaS.  isi(2)  .Ni.o')  GO  TO  '^'S 
L8  ■ L8  + I 

WRITE  (lh384)  T(l).  P.  DP.  RATE 
384  FORMAT  (4E13.6) 

675  CONTINUE 
T22  - T(l) 

L9  ■ L9  + I • 

IF  (T(l).GE.TSTOP)  N - 0 
IF  (N,EQ.0)  WRITE  (11,384)  XY2.XY2,XY2.XYZ 
RETURN 

100  FORMAT  (IHI/SX,5A4/7X,'TIME',6X,*PRESS',6X  ,'DP/Dr,eX, 

iilt‘,bX,'UI  KR'«4X,  bUKh  t-k 'Uibb  tikli  ,bl^, 

2 ' PDOT/PMAX' /7X. * MSEC* .SX, * KPS  I A' , 4X. ' MPS IXSEC' ,5X, 

3 'IN/SEC',27X,' INCH', 6X,* MSEC- 1*/) 

101  FORMAT  (2(4X.F7.3).4X,F8.3.4X.F7.3,3(4X,F7.4).4X,F8.3) 
END 
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0301 

0002 

0003 

0004 

0005 

0006 
000? 

0000 

0009 
0011 
0012 

0013 

0014 

0015 

0016 
0017 

0010  10 

0019 
0021 
0022 

0023 

0024 

0025 

0026 
0^27 

0020 

0029  300 

0030 

0031 


SUBROUTINE  FORMKDB,  I,0D1.PD1,UID,U0D,  V0,GL.S.UCRIT) 
DATA  P 1/3. 141593/, P 14/. 785399/ 

U • DB>t<2.0 

DDUt  - (UOD  - UID)/4. 

QD  ■ ODl  + DDL) 

PD  • PDl  - DDU 
PD  - AMAXl  (PD,  0.0) 

OD  - AMAXl  (OD,  0.0) 

IF  (I.NE.0)  GO  TO  10 
E • PI4)*<(0Diiot!2  - PDiM«2) 

E • AMAXl  (E,  0.0) 

V0  - E*GL 

S ■ 2.*E  GLn<PI>x(OD  + PD) 

I • I + I 

UCRIT  « (OD  - PD)/2. 

RETURN 

UCRIT  - (OD  - PD)/2. 

IF  (U.GT. UCRIT)  GO  TO  300 

OOD  - OD  - U 

ODD  - AMAXl  (OOD,  0.0) 

PPD  • PD  + U 
GGL  - 6L  - U 

E - PI4'K(Q0D*>i<2  - PPD**2) 

E • AMAXl  (E.0.0) 

S - 2.*E  + GGL*PI#(00D  + PPD) 

RETURN 
S - 0.0 
RETURN 
END 
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0002 

0003 

0004 
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0006 
0000 

0009 

0010 
0011 
0012 

0013 

0014 

0015 

0016 
0017 
0010 

0019 

0020 
0021 
0022 

0023 

0024 

0025 

0026 
0029 

0029 

0030 

0031 

0032 

0034 

0035 

0036 

0037 

0038 
@039 

0040 

0041 
0S43 

00.14 

0046 

0947 

0048 

0850 

0051 

0052 

0053 

0054 


FORTRAN  IV 

SUBROUTINE  FQRM7(DB. I,D,PD,UID,UOD,V0,GL,S.CRITU) 

LOGICAL  RYPE 

DATA  PI,PI3,PI4/3.14 159265359 . 1 . 047 1 9755 1 1 97, . 785398 1 633975/ 
DATA  RT/l. 7320508076/ 

U ■ 2.>*<DB 
IFd.GT.O)  GOTO  10 
I - 1 

E -P 1 4*  ( D'M‘2-7 . ♦PDiM<2 ) 

VO  ■ 

V-V8 

S0»2 . *E+P  l>t<(D+7 . *PD)  x-GL 
S«S0 

Uy-PD+UID 

E- . 5h<(  Dxok2-PD>k»2+4  . i«ULJhok2-2  . )rtJUx<Dx<RT)  / ( D+PD-UU*RT) 

CRITU“AMINl  ( .5*(D-PD)  ,AMAX1(2.*IJUM?T-PD.E)  ,GL) 
RYPE-D-UlD.LT.UUx-RT 
2-0. 

RETURN 
10  E-0, 

UMIN  - AMINl  (yiD.UOD) 

RYPE-D~yiD.LT.yUx<RT 

S-0. 

U-AMAXl(U.8.) 

GRL-GL'U 

IF(U.GT.CRITU)  GOTO  300 

prfd-pd+u 

PRFD2-PRFD**2 

OD-D-U 

002-OD>Kit<2 

IF(U.GT.miN)  GOTO  100 
E-Pl4x<(0D2-7.ii<PRFD2) 

S -2 . x€+P  I x;  (QD+7  . x4>RF0)  *GRL 
GO  TO  309 
100  LSJ-PD+yiD 
yu2-yy**2 

THETA-2 . *ACOS (AMIN  I (LU/PRFD. I . ) ) 
ALPHA-ACOS(RMtNl((.25x<(OD2-PRFD2)+Uy2)/(yy*OD).  1.)) 

IF(ALPHA.GE. .5»PI3)  GOTO  218 

BETA-ACOS(AMAXi  ( ( .2S>»*(PRFD2-0&2)+yLi2) /<yy*PRFD) 1 . ) ) 

I -.SxcTHETA  - P13 
IF(RVPE)  GOTO  250 

E-3 . I N (ALPHA)  H . S*( 0D2»( . 5x«P  1 3-ALPHA)  -La.J2«RT-PR?D2*<0ETA 

I x-.5*S  IN  (THETA))) 

S -2 . *£ + ( 6 . »SE  TA*PRF  ( P I -6 . *ALPH0 ) xcOD ) xGRL 
210  1F(THETA.GE,PI3)  GOTO  300 

E 1 - 1 . 5*(Uy2*8T- 1 . 5*PRFD2«(S IH(TWETA) +P  13-THETA) ) 

E-E+El 

S-S+2 . #£ I +9 . *PRFD*(P 1 3-THETA >*GSL 
GOTO  300 

250  £-3.»UJ*t)Dx^SINCAiPHA)+l.S*(0D2*(.3«?I3 

1 -ALPHA)  -PRFD2*(SIN(THETA) 

2 +eETA*.5*(Pl-TKETA))) 
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0055  S-2 . *e+6 .*(0D*( , 5*P I3-ALPHft) +PRFD*(BETA+. 5*(P ! -THETA) ) ) ♦GRL 

eass  300  v-E^rtiRL 

0057  2-1.-V/V0 

0058  RETURN 

0053  END 
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0031 
0932 
0033 
0834 
0035 
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0048 

0041 

0042 

0043 

0044 
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0046 
9040 

0049 

0050 
0851 

0053 
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0055 
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SUBROUTINE  F0RM19(DB, I.D>PD>UID,UDD, V0.GL.S,CRITU) 

LOGICAL  RYPE 

DATA  P LP13.P  14/3. 14159265359,  l.04719?S5n97..?83398l6339?S/ 
DATA  RT/l.?32050a076/ 

U ■ 2.H<D0 
IFd.GT.a)  GOTO  10 
I - I 

E-P  I4*(:D>toti2- 19.ikPD>k>k2) 

V0-E>«3L 

V“V0 

sa -2 . >i€+P  I ♦ (D+ 1 9 . ♦PD)  >t<GL 
SoSO 

UU^'RO+UID 

E“(D4fD)>to«2 

UU2-UUW2 

E-.5><<(D-PD-ULJ;HSQRT(12.»(E-I6.!»!ljy2)/'(E-  ,2.>WJU2))) 

CR I TU -AMI  N U AflAX  1 ( 2 . >t<UU4?T-PD,  E ) , GL ) 
RyPE-D-UID.LT.UU*SQRT(13.) 

IF (RYPE)  GOTO  250 
2-0. 

RETURN 
10  E-0, 

UMIN-AMINKUID.UOD) 

S-O. 

U-AMAXKU.O.) 

CRL-GL-U 

rF(U.GT.CRlTU)  GOTO  30S 

PRFO-PD+U 

PRFD2-PRFD»=*t2 

OD-D-U 

0D2-0D=*‘^f2 

IF(U.GT.UMIH)  GOTO  100 
E-PI4*{0D2“19.*PRFD2) 

S - 2 . ♦£ ♦P I *( OD + I 9 . ♦PRFD ) oCRL 
GOTO  300 
100  LU-PDH)ID 
UU2-UU«*2 

■n !E TA - 2 . !»«COS  ( AMI M 1 ( LU/PRF D * I . ) ) 

ALPHA  I -ACOSC AMIN  1 ( ( . 12S«<(0D3-PRFD2)+2. ♦UU2)  /(IJJ*0D) . 1 , )) 
ALPHA2-AC0S(AM1H  I ( ( ,25»(OB2-PaFD2)+3.*LIJ2)/(LSJ«OD'*«?T) , 1. ) ) 
ALPHA-ALPHA I ♦ALPHA2 
IFCALPMR.GE. .5»P13)  GOTO  210 

BETA  I -ACOS  ( AMAXl  ( ( . l25*(PRFD2-0D2)+2.  »L!U2)  /(LU-PRFD)  ,-!.>) 

I -.5»THETA“Pi3 

BETA2  - ACOS  ( AMAX I ( ( , 25*  (PRFD2-OD2 ) +3 . *Uy2 ) / (liJtF'RF  B*8T) . - 1 . ) ) 

I -.5«(THErA+P[) 

BETA-BETA  1+45ETA2 
IF (RYPE)  GOTO  259 

E -3 . *0D*5JU*  { 2 . liS  1 N ( ALPHA  J ) ♦RT*S  I N ( ALPHA2 ) > -6 . ‘fyy2*R  T*  1 . 5* 

I (0D2*(.S'^I3-ALPHn)-PRFD2-»'(SlN(TKETA)+{3ETA)) 
S“2.»€+6.*(0D*(  .5*PI3-ALPHA)'K»RFD48ETA)>*GfiL 
210  IF(THETA.GE.PI3)  GOTO  303 
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0057  Et-6.*(lJLJ2*RT-1.5*PRFD2>«(Sm(THETA)+PI3--THETfl>) 

0058  E-E-fEl 

0059  S-S+2 . *E 1+36. »PRFD*(P I3-'mETfl) *GRL 

0060  GOTO  300 

0061  250  TYPE  251 
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0062  251  FORMfiTC  ' GRAIN  GEONETliY  YIELP*  NO  OUTER  SLIVERS.') 

0063  STOP 

0064  300  V-E>t«RL 

0065  2-1.-V/V0 

0066  RETURN 

0067  END 


RT-ll 
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0805 
0886 
000? 
0809 
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6014 

0015 
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FOSTRflN  IV  Vei-1^  SOURCE  LISTING  PAGE  001 

SUBROUTINE  FORMSP  <D8, I,ODl.PDl,yiD,UOD,V0>GL,S,UCRIT) 
BftTA  PI4^-.r85399/.PI/3. 141393/ 

U - DBX‘2.0 

DDU  • (UPD  - UID)/2. 

OD  - ODl  + DDU 
OD  - AMAXl  (OD,  0.00 
IF  (I.NE.0)  GO  TO  10 
I - I + I 
R - UD/2,0 
S “ 4.*P!»Ri*3!<2 
va  - 4.*PI/3.»Rkc)K3 
UCRIT  - OD 
RETURN 

18  UCRIT  " OD 

IF  (U.GT.UCRrn  GO  TO  300 
R - (OD  ~ U)/2. 

S ■ 4.>i‘PI>«R>w<2 
■RETURN 
300  S - 0.0 
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RT-U 

0901 

0002 

0033 

0004 

0005 
0036 
0807 
0038 
0009 
60  U 
00J2 
0013 
a6l<s 

0015 

0016 
8317 
6018 
0020 
0021 
0022 
0023 
3024 

0025 

0026 
0027 
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SUBROUTINE  FORMrV  (.JS,  I.ODhPD,UID,UCD,V0.GLUS,UCRIT) 
DATA  P 14/. 7a5399/.PI/3. 141593/ 

U " D0W2.8 

DDU  - (UOD  “ UID)/2. 

OD  » QDl  + DDU 
GL  - GLl  + DDU 
OD  - AMAXl  (OD,  8,0) 

GL  - AMAXl  (GL,  8.0) 

IF  (I.NE.0  ) GO  TO  10 
1 - I + 1 
E « PI4!»'-0D#*2 
va  - E>«GL 

S - 2. IKE  + PIiKGLiKOD 
UCRIT  • AMINl  (QD,  GL) 

RETURN 

10  UCRIT  - AMINl  (QD,  GL) 

IF  (U.GT. UCRIT)  GO  TO  303 
OOD  - OD  - U 
GGL  ■ GL  ~ U 
E “ Pt4!K00D'M<2 
S • 2,*E  + PI#8GLiOOD 
RETURN 
300  S « 0.0 
RETURN 
END 
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0001 

0Q02 

0003 

0004 
0066 
000? 

6088  10 
0008 
0010 
0012 
0813 
8014 


FUNCTION  flCOS  (ft) 

DftTft  PI2/1. 570796/,  PI/3.141593/ 
B » SORT  (1.  - ft**2) 

IF  (A.NE.0.0)  GO  TO  10 
ftCQS  - PI2 
RETURN 
C • B/A 
n - ATAN  (C) 

IF  (A.LT.e.0)  D - D + PI 

ACQS  - D 

RETURN 

END 
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SUBROUTINE  UALTEM  (F>DTINE,DUM1.DUM2) 

COMMON  RUN  ID (5) , RTITLE (18), DATE (3) . OPERN (S) , PROPN ( IB) , 

1 PSORCEC 18)  ,PRl.OT(  18)  ,PREM(  18)  ,TIGNR{  18)  ,GAGE<  18) , 

2 SPACE  1(36), CHGUT, UID, PTEMP. UGHTI , PCORR, BVOL, BTEMP, 

3 SCAL,CALrnN,STIME,FFID,BL,OD.PD,UOD,Fl,DUM,XMUXl,ETl,GAMU 

4 SPACE2<4), I3K1,KTQT,M2,MY,PMAX,TMAX,NPMA,PPM,DPMAX,IHL,H, 

5 RHaC,RHa,Tl0,T90,TiaS0>P9(S),DP9(5),NPTH,PX(20), 

6 FX(20),ETX(20),XMUX(20),GAMX(20), IDFF,NPA,UEBX(20),ABX(20), 
? ISK, ISK2,MEM2,XMI,RP.T2:ER0,V0L,CP,CST,TM,PIGN,TIGN,PTHE0, 

0 CaNl,C0N9,M£M3,MEM7,SFA2.L9.P,DP,TST0P, 19, 

9 L8, CONU, TBU, FRAC, 2, SRAT, XTBU, RATE , RL 
COMMON  /KEVIN/  TE(35),  XT(35) 

DIMENSION  TMA(2),  TTA(2),  DELX(35),  DTt35) 

D ATA  I SUS^'TI/,  ALF A/ . 0 109 1 /,  XK/ . 000666?/,  DELX0/ . 0009/,  FAC/0 . 3/ 

TTA(l)  - DUMl 

TTA(2)  - DUM2 

TMA(l)  - 0.0 

TSUM  - 0.0 

DTIME  ■ DTINE/1000. 

TMA(2)  • DTIME 
IF  (ISUS.EQ.G)  GO  TO  86 
IF  (DELTCN  - DTIME)  90,90,91 

91  DELTIM  » DTIME 
GO  TO  92 

90  DELTIM  ■ DELTCN 

92  TSUM  - TSUM  + DELTIM 

IF  (TSUM  - DTIME)  93,93,95 
95  TSUM  • TSUM  - DELTIM 
DELTIM  - DTIME  - TSUM 
TSUM  - TSUM  + DELTIM 

93  CONTINUE 
MEM?  - I 

IF  ((TMA(2)-TMA(1)) .EQ.0.0)  GO  TO  1942 
CALL  FIND  I (TSUM,TG.Tt«,TTA,2,rEM?) 

GO  TO  1943 

1942  TG  - TTA(l) 

GO  TO  1943 

1943  QDOT  • F*(TG  - XT(1)) 

ODOT  - QDOT  + SPACE2(3)*(TG*>wl  - XT(l))Wt4) 

T2ER0  - QDOT/XK#2.>kDELX0  + XT (2) 

DT(1)  • ALFA/(DELX0x«k2)>k(T2ERO  - 2.>tO<T(l)  XT(2)) 

DO  42  I * 2,31 

CONI?  ■ C0NA/(DElXlI)>K>f2) 

42  DT(I)  - C0N17H<(XT(l-l)  - CONB^XTd)  + CONC*XT<I+U) 

DO  S I • 1,31 

5 XTd)  • XTd)  + DELTlM+DTd) 

XT<32)  - XT(30) 

IF  (TSUM  - DTl?E)  92,74,74 
74  DELTIM  - DTIME 
RETURN 
86  ISU6  « 1 

DO  21  I • 1.35 
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6046  >CT(n  - 296. 

0047  21  TE(I)  - 298. 

0048  TYPE  777,  ALFA 

0049  777  FORMAT  (5X.'80MB  THERMAL  DIFFUSIVITY',F12.5) 

0058  ACCEPT  778,  TEtf> 
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0051 

0052 

0054 

0055 

0056 
IIQ57 

nqsa 

0059 

i-'»i60 

oOei 

0062 

0063 

0064 

0065 


778  FORMAT  (E12.0) 

IF  (TEMP. NE. 0.0)  ALFA  • TEMP 
CONA  - ALFA/d.  + FAC/2.) 
CONB  - (2.  + FAC) /(I.  + FAC) 
CONC  - l./(l.  FAC) 

J - 1 

DUEB  - DELX0 
DELX(J)  - DELXa 
DO  22  J ■ 2,34 
DELX(J)  • DUES 
22  DUEB  - DUEBixCl.  + FAC) 

DELTCN  • 0.35>i‘DELX0iW<2/ALFA 

RETURN 

END 
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0001 

SUBROUTINE  RCALC 

C 

SUBROUTINE  TO  DO  LEAST  SQUARE?  FIT  ON  BURN  RATE  DATA 

0002 

COMMON  RUNID(5),RTITLE(18).DaTE(3),0PERN(5),PP0PN(18). 

1 PSQRCE( 18) ,PRLDT( 18) ,PREM( 13) .TIGNR( 18) .GAGE( 10) , 

2 SPACE  1 (36) > CHGUT,UID, PTEMP. UGHTI , PCORR, BVOL. BTEMP, 

3 SCAL,CALCON.STIME,FPID,GL.OD,PD.UOD.Fl,DUM.XrilJXUETl.GAMU 

4 SPACE2(4).  ISKl,KT0T,(12,M'f'.PMAX,TMAX,NPMA,PPM,DPMAX,  IHL.H, 

5 RHOC,RHO,T10.T90,T1090.P9(5),DP9(5),NPTH,PX(20), 

6 FX(20).ETX(20).XMUX(20).GAMX(28),IDFF.NPA-LEBX(2O),A8X(20), 

7 ISK, ISK2.MEM2>XMI,RP,T2ERO,VOL*CP,CST.TM,PIGN.TIGN.PTHEO, 

8 CON1,CON9,MEM3>MEM7,SFAC,L9,P1>DPUTSTOP.  19. 

9 LS.CQNU.TBy.FRAC.Z.SRAT.XTBU 

0003 

DIMENSION  T(400),P(400),DP(400).R(4a0) 

0004 

COMMON  /BLAH/  PLQ.  PHI.  ACO.  XNC.RSQ.  Rl^ 

6005 

REUIND  11 

0006 

DO  669  I - 1.L8 

0007 

READ  (11.384)  T( I) .P( I) .DP( I) .R( I) 

6000 

384 

FORMAT  (4E13.6) 

0003 

669 

CONTINUE 

6010 

PLL  - PLO>*<PMAX 

0B11 

PUL  - PHIi*<PMAX 

0012 

K - 0 

0013 

SUMl  - 0.0 

0014 

SUM2  - 0.0 

0015 

SUM3  - 8.0 

0016 

SUM4  -0.0 

0017 

SUM5  • 0.0 

0018 

DO  735  I • l.LB 

0019 

IF  (P(I).LT.PLL)  GO  TO  735 

0021 

IF  (P(I).GT.PUL)  GO  TO  735 

0023 

K - K + 1 

0024 

OUMl  • ALOG  (P(I)*1000.) 

0025 

DUM2  » ALOG  (AMAXl (R( I ) ,8.001 )) 

0026 

SUMl  - SUHl  + DUMl*DU«2 

0027 

sure  • suti2  + dumi 

0028 

SUM3  - SUM3  + DUM2 

0023 

SUM4  - SUM4  ♦ DUMl>toK? 

0030 

SUMS  » SUMS  -s-  DUM2**2 

0031 

735 

CONTINUE 

0032 

XNC  - (SUMl  - SUM2*SUM3/K)/(SUM4  - (SUtt2iMi2)/K) 

8033 

ACO  - EXP  (SUM3/K  - XHC>t<SUM2/K) 

0034 

RSQ  - ((SUMl  - SUM2<SUM3/K)>m<2) '((SUM4  - (SUMaiM^l/K)* 

— ' 

1 (SUMS  - -(SUM3**2)/K))  ■ 

0035 

SUM-  8.0 

8036 

K • 0 

6037 

DO  733  I • I.L0 

■ 0039 

IF  (P(I).LT.PLL)  GO  TO  733 

0048 

IF  (P(!).GT,PUL)  GO  TO  733 

■ 9042 

K • K + 1 

0043 

DUMI  • ACO*(P(n*10O0.)**XNC 

0044 

DUMI  • (R(I)  - DUMI )/AMAXl(R( I). 0.001) 

0045 

SUM  • SUM  + DUMl>M(2 
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0046 

0047 

0048 
0043 


733  CONTINUE 

RMS  - SORT  (SUfV<K  - 2))»100.0 

REIURN 

END 
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C VERSATRAN  PLOTTER  SUBROUTINE  FOR  CBRED 

0001  SUBROUTINE  VPLOT 

0002  COMMON  RUNID(5),RTITLE(1B),DATE(3),OPERN(5),PROPN(I0), 

1 PSQRCE(l0).PRLnT(10),PREM(18),TIGNR(I8),GAGE(lB), 

2 SPACE  I (36) . CHGUT, UID, PTEMP,UGHTI > PCORR, BVOL, BTEMP, 

3 SCAL,CALCON>STIME,FPID,GL,OD,PD,UOD,Fl.DUM>XMUXl,ETl.GAMl. 

4 SPACE2 (4) . ISK 1 , KTOT, M2, MY, PMAX, TMAX, NPMA . PPM. DPMAX, IHL, H, 

5 RHQC,RHO,T10,T9e,T1090,P9(5),DP9(5).NPTH,PX(20), 

6 FX(2a),ETX(2B),XMUX(20),GAMX(20),IDFF.NPA,UEBX(20).ABX(20), 

7 ISK,  1SK2,MEM2.XMI.RP,TZER0,V0L.CP,CST,TM.PIGN.TIGN,PTHEQ, 

8 C0Nl>C0N9,MEM3,MEM7.SFAC,L9,Pl,DPl,TST0P, 19, 

9 L8,C0NU,TBU,FRAC,2,SRAT,XTBU 

0003  DIMF.NSION  T(408) ,P(400)  ,DP(400) ,8(400) 

0004  COMMON  /BLAH/’  PLO,  PHI.  ACO,  XNC,  RSQ,  RMS 

0005  DIMENSION  ABC (17) 

0006  DIMENSION  BCB(15) 

0007  DIMENSION  A(ll) 

0008  . DIMENSION  T9(2),  PI  1(2) 

0009  DIMENSION  Xl(9),YI(13) 

0010  DIMENSION  TIT(16) 

0011  DATA  TIT/* P/PM', 'AX  ','PDOT','  '.'RUN  ','IDi  ','RUN  ', 

1 'TITL'.'E:  '.'PROP*,'  TYP'.'Ei  '.'GRAI'.'N  TY'.'PE:  '. 

2 'AT  '/.  S/9999./ 

0012  DATA  8/1.0,2.0,4.0.6.0,0.0,10.0,20.0.40.0.60.0,80.0,108.0/ 

0013  DATA  Xl/'TIME'.'-MSE'.'C  '.4HP/PM.'AX  ','1.0 

1 *10.0'. '100. '/,Y1/'PRES'.'S-KP','SIA  ','PDOT','  MPS', 

2 4HI/SE.'C  ','0.1  '.'UO  '.' 10.0*. 'RATE*. 4H-IN/.'SEC  '/ 

0014  DATA  AQC/'THE  ', 'CONS' . 'TANT' . 'S  IN'.'  THE'.'  EQU'.'ATIO*. 

1 'N  R A*'.'P**N'.'  ARE', 'A!  ','N:  ','FOR  '.'P/PM', 

2 'AX  ','  TO'/ 

0015  DATA  BCD/'COEF'.'FICI'.'ENT  '.'OF  D' , 'ETER' . 'MINA' , 

1 'TION'.'  : '.'PER  ','CENT*,*  ROO' , 'T  ME' . 'AN  E'.'RROR'.'  $ '/ 

0016  REWIND  11 

0817  DO  683  I - l.LB 

0016  READ  (11.304)  T( l).P( I) .DP( I) .R( 1) 

0819  3B4  FORMAT  (4E13.6) 

0020  603  CONTINUE 

8821  ENDFILE  1 1 

0822  CALL  MODE  (0,2.0, 1.0, -I .) 

0023  CALL  MODE  (2,0.0,-0.75,0.873) 

0024  CALL  MODE  (3, 6. 0,-0. 75. 3. 8) 

0025  CALL  MODE  (7.S.0.5.0.S) 

0026  TYPE  1976 

0027  1976  FORMAT  (5X.'D0  YOU  WANT  TO  SUPPRESS  THE  IGNITION'/ 

1 5X,'TRC  LAG  ON  THE  PRESSURE  AND  DP/DT  PLOT?'/ 

2 SX.'TVPE  "S"  TO  SUPPRESS') 

0020  ACCEPT  1977.  STS 

0029  1977  FORMAT  (A4) 

0030  IF  (STS.NE.STR)  GO  TO  1978 

0032  DATA  STR/'S'/ 

0033  T9(l)  - lO.*rFIX(TIGN/ie.) 

0834  DUMA  • 0.0 
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0035  GO  TO  1979 

0036  1970  T9(l)  - 0.0 

0037  DUMfl  - 5.0 

0038  1979  T9 (2)  - TMAX 

0039  PlUl)  • 0.0 
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Pll(2)  • PMfiX 

0041 

CALL 

SCAN 

(T3,  Pll,  -2.  440) 

0042 

CALL 

MODE 

(-8,DUMUDUM2,DUM3) 

0043 

DUK2 

- AMAXl(DUM2,  DUMA) 

0044 

CALL 

MODE 

(0.  DUMU  DUM2>  S) 

0045 

DUM4 

- DUM2 

0046 

DUMB 

- DUMl 

0047 

CALL 

MODE 

(-9.  DUMl,DUM2>DUM3) 

0048 

DUf12 

- AMAXl(DUM2,5.0) 

0049 

CALL 

MODE 

(9.0.0.DUM2,S) 

0050 

CALL 

AXES 

(9.2.Xl(l),n.2,Yl(l)) 

0051 

CALL 

DRAU 

(T.  P,  L8.  440) 

0052 

CALL 

SCAN 

(T,  DP.  -LB.  440) 

0053 

CALL 

MODE 

(8.DUM8.DUM4.S) 

0054 

CALL 

MODE 

(-2.DUMl,DUM2.DUM3) 

0055 

DUM5 

- DU^O  - 0.5 

0056 

CALL 

MODE 

(2.  S.S.DUM5) 

005? 

CALL 

raDE 

(9.  0.0.  S,  S) 

0056 

CALL 

AXES 

(-4.0,X1(4). 13.3.Yl(4)) 

0059 

CALL 

MODE 

(2.S,S.DUM3) 

0060 

CALL 

DRAU 

(T.  DP.  LB.  440) 

0061 

K - 0 

0062 

CALL 

MODE 

(4, .073. .055, S) 

0063 

CALL 

MODE 

(6.S. .080.S) 

0064 

CALL 

MODE 

(3.S,S.-.2) 

0065 

CALL 

NOTE 

(0.0.3.0.T1T(5).7) 

0066 

CALL 

NOTE 

(!.0.3.0.RUNiD(l).20) 

0067 

CALL 

NOTE 

(0.0.2.0.TIT(?).l0) 

0066 

CALL 

NOTE 

(1.0.2.8.RTITLE(1),60) 

0069 

CALL 

NOTE 

(8.0.2.6.T1T(10), 10) 

0970 

CALL 

NOTE 

<l.B.2.6.PRQPN(l).6e) 

0071 

CALL 

NOTE 

(0.0.2.4.TiT(l3).  11) 

0072 

CALL 

NOTE 

(1.0.2.4.FFID.4) 

0073 

IF  (K.EQ.0>  GO  TO  10 

0075 

CALL 

NOTE 

(0.0. 2. 2. TITO). 4) 

0876 

CALL 

NOTE 

(2.0.2.2.TIT(16).2) 

0077 

CALL 

NOTE 

(4,0,2.2.TJT(l).6) 

0079 

Y • 2 

[.2  - 

0.2 

0079 

DO  61 

1 1 •* 

1.4 

0080 

CALL 

NOTE 

(0.0.Y.DP9( I). 1003) 

0081 

CALL 

NOTE 

(4.0.Y.P9(l).  1003) 

6002 

61 

Y - Y - 0, 

2 

0003 

10 

CALL 

MODE 

(4,.  U.067.S) 

0084 

CALL 

MODE 

(6.S. .l.S) 

0085 

CALL 

MODE 

(3.S,S,3.0) 

0006 

CALL 

DRAU 

(0..0.. 1.9080) 

8897 

pmx 

- 0.0 

0088 

DO  65 

i 1 • 

1,L8 

0083 

65 

PtIAX 

- AMAXl(PMAX.P(l)) 

6030 

DO  69 

1 I -1 

.L8 

0091 

69 

T(l) 

- P{l)/PMAX 

8892 

CALL 

MODE 

(Q.  0.0.  0.2.  0.0) 
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0093  CfiLL  AXES  (6.2,Xl<3), 13.3-Yl<4)) 

0094  CALL  DRftU  (T,  DP*  L0*  440) 

0095  CALL  FORM  (5*  1.0*  S.  1.0) 

0096  K ■ t 

0097  CALL  MODE  (4* .073* .055*3) 
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0098 

0099 

0100 
0101 
0102 
0183 

0104 

0105 
0186 

0107 

0108 
0110 
0111 
0112 

0113 

0114 

Oils 

0116 

0117 

one 

0119 

0120 
0121 
0172 
0123 

0125 

0126 

0127 

0128 

0129 

0130 

€131 

0132 

6133 

613S 

0136 

0137 

0138 

0139 

0143 

0141 

0142 

0144 
0146 
814? 
0148 

0143 

0150 

0151 

0153 

0154 


CALL  MODE  (6,S. .080.9) 

CALL  MODE  (3.S.S.-.2) 

CALL  NOTE  (0.0,3.0.T1T(5).7) 

CALL  NOTE  ( 1 .0.3.0.RUNID ( 1) .20) 

CALL  NOTE  (0.0.2.8.TIT(7) . 10) 

CALL  NOTE  (1.0,2.8.RTITLE(l)-6a) 

CALL  NOTE  (0.0.2.6.TIT(10), 10) 

CALL  NOTE  ( 1 .0.2.6.PROPN( 1) .68) 

CALL.  NOTE  (0.0.2.4,TIT(13).  11) 

CALL  NOTE  ( 1 ,0.2.4.FFID.4) 

IF  (K.EQ.0)  GO  TO  11 
CALL  NOTE  (0.0.2.2.TIT(3).4) 

CALL  NOTE  (2.0.2.2.TIT( 16).2) 

CALL  NOTE  (4.0.2.2.T1T(1).6) 

Y - 2.2  - 0.2 
DO  91  I - 1.4 

CALL  NOTE  (0.0.Y.DP9(1). 1093) 

CALL  NOTE  L4.0.Y.P9(I).1803) 

81  Y « Y - 0.2 

11  CALL  MODE  <4. . 1. .867.5) 

CALL  MODE  (6.S..1.S) 

CALL  MODE  (3.3. 5. 3. 9) 

CALL  DRAU  (0. .0. . 1.3000) 

PFAC  - 1.0 

IF  (PMAX.LE.10.0)  PFAC  - 10,0 

PMIN  - 1.0 

PM  IN  - PMIN/1>FAC 

IGQ  • 0 

PSTO  • PMAX 

TVPE  S'^9 

529  FORMAT** (5X.‘ DO  YOU  UANT  ALL  THE  RATE  CURVE?'-' 

1 5X. ‘ENTER  I FOR  YES.  PLEASE') 

ACCEPT  530.  IGO 

330  FORMAT  (11) 

IF  (IGO.EO.l)  GO  TO  533 
PMIN  • PLO»PrWX 
PSTO  - PH1*P(1AX 

333  CONTINUE 

RMAX  - 10.0 

RMIH  • 0. 1 

J97  - 0 

00  69  I - I.L9 

IF  (Pd)  .LT. PMIN)  GO  TO  60 

IF  (Pd), GT. PSTO)  GO  TO  68 

J97  - J37  + 1 

P(J97)  • ALQG19  (Af^Xl (P(  1 ) .PMIN) ) 

R(J97)  • ALOGIO  (AMAXUAHIHURd) . 10.6).  0.1)) 
68  CONTINUE 
DUMl  • 0.0 

IF  (PFAC. GT. 1.0)  DUMl  • -1.9 
CALL  lODE  (8. DUMl,  0,4.  S) 

CALL  MODE  (9.  -1.0.  0.4.  S) 
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0155 

0156 
015? 
0158 
0160 


DO  66  1 ■ 1.10 

66  AH)  • 2.S#ALOG10(A(I+l)/A(I)) 
CALL  FORH  ( 1010.A. 1610.A) 

IF  (PFAC.EQ. 1.0)  GO  TO  997 
CALL  NOTE  (-0.  l.-0.25,YU0).3) 
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0161 

0162 

0163 

0164 

0165 

0166 
016? 
0166 
0169 
01?0 
R171 

0172 

0173 

0174 
0173 

0176 

0177 

0178 

0179 
0160 
0181 
0182 

0183 

0184 

0185 

0186 

0187 

0188 

0189 

0190 

0191 

0192 

0193 
6194 

0195 

0196 

0197 

0198 
0159 
0200 
0201 
0262 

0203 

0204 

0205 

0206 
0207 

0209 

0210 
0211 
0212 
0213 


CftLL  NOTE  (2.3*-0.25,Yl(3),3) 
CftLL  NOTE  (4.7,-0.25,YlC10),4) 

GO  TO  998 

997  CftLL  NOTE  l.-0.25,Xl(7),3) 
CftLL  NOTE  (2.3>-0.25,  XI (8),  4) 
CftLL  NOTE  (4.7.-0.25,  XI (9),  4) 

998  CONTINUE 

CftLL  flQDE  (4,  S,  S>  90.0) 

CftLL  NOTE  (-.15>-.10>  Yl(8),  3) 
CftLL  NOTE  (-.15,2.40,  Yl(9),  3) 
CftLL  NOTE  (-.15,4.8,  Yl(10),  4) 
CALL  HQDE  (-4,DUMl,DU!12,DUt13) 
CftLL  MODE  (4,  .t6,.10,S) 

CftLL  MODE  (-6,DUM4,DUMS,DUIiS) 
CALL  MODE  (6, S. . 16, .213) 

CALL  NOTE  (-.35, 1 .75.Y1 ( 1 1) , 1 1) 
CftLL  MOPE  (4,  S,  S,  0.0) 

CALL  NOTE  (1.90,-. 51,  Yl(l),ll) 
''ALL  MODE  (4,DUM1,DUM2,S) 

CALL  MODE  (6.S,DUM5.DUM6) 

CALL  DRAU  (P.  R,  J97,  442) 

K * 8 

CALL  MODE  (4, .073, .055.S) 

CALL  MODE  (6,S,.080,S) 

CALL  MODE  (3,S,S,-.2) 

CftLL  NOTE  (0.0,3.0.TIT(5),7) 

CftLL  NOFE  (1.0.3.0,RUNID(l),20) 
CALL  NOTE  (0.0,2. 8, TIT(7), 10) 
CftLL  NOTE  (1.0,2.0,RTITLE(1),68) 
CftLL  NOTE  (0.0.2.6,TIT(1O), 10) 
CftLL  NOTE  (l.0,2.6,PROPN(n,68) 
CALL  NOTE  (0.0,2.4,TIT( 13) . 11) 
CftLL  NOTE  (1.0,2.4,FFID,4) 

CftLL  NOTE  (0.0,2.2,ABC(1),44) 
CftLL  NOTE  (0.0,2.0,ftBC(12),3) 
CftLL  NOTE  (1. 0,2. 0,ACO, 1006) 

CALL  NOTE  (0 . 8. 1 .8.ftBC( 13) . 3) 
CALL  NOTE  ( 1 .8, I .8.XNC, 1003) 

CALL  NOTE  (0.0, 1 .6,ftBC( 14) , 11) 
CftLL  NOTE  (l.0,l.6.PLO, 1003) 

CftLL  N0TE(1.4,l.6.flBC(l7),4) 

CALL  NOTE  <1.72,  l.G,PHI, 1003) 
CALL  NOTE  (0.0. l.4,BCD< 1) .31) 
CALL  NOTE  (2.43, 1.4, RSiJ.  1004) 
CftLL  NOTE  (0.0.1.2,BCD(9).27) 
CALL  NOTE  (2.48,1.2,8115,1004) 

IF  (K.EL.O)  GO  TO  12 
CftLL  NOTE  (0.8.2.2,TiT(3).4) 

CALL  NOTE  (2.a.2.2,TIT( 16) ,2) 
CALL  NOTE  <4.0,2.2, TIT( 1) ,6) 

Y - 2.2  - 0.2 
DO  SI  I - 1,4 
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0214  CALL  NOTE  (0.0^Y,DP9(n,  1003) 

0215  CfiLL  NOTE  (4.0,Y.P9(I). 1003) 

0216  91  Y - Y - 0.2 

0217  12  CfiLL  MODE  (4, -U. 067,8) 

0210  CfiLL  MODE  (P.S,.U3) 


t 
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0219 

CALL 

MODE 

0220 

CALL 

DRAU 

0221 

CALL 

DRAU 

0222 

RETURN 

0223 

END 

(3.S,S,3.B) 

(e./B., 1.9000) 

(6.0.  e.0.  0.0.  9999) 
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0001  SUBROUTINE  FINDTP  (ARG.  ANS,  NVAR,  NARG,  NANS,  IDEV,  NPTS,I,M3) 

C SUBROUTINE  FIND.  TABLE  LOOK-UP  UITH  LINEAR  INTERPOLATION. 

C FINDTP  IS  A DIRECT  ACCESS  READ  LOOKUP  MODIFICATION  OF 
C FIND.  INSTEAD  OF  AN  ARRAY,  PASSING  THE  ARGUMENTS,  IT 

C EXPECTS  A FILE  WHICH  HAS  BEEN  FILLED  UITH  NVAR  VARIABLES, 

C TO  THE  EXTENT  OF  NPTS  RECORDS.  THE  INDEPENDENT 
C VARIABLE  MAY  OCCUPY  ANY  OF  THE  LOCATIONS  l-NVAR  AND 
C NARG  GIVES  ITS  LOCATION.  LIKEWISE  FOR  THE  DEPENDENT 

C VARIABLE  WHOSE  LOCATION  IS  GIVEN  BY  NANS.  IDEV  GIVES 

C THE  DEVICE  NUMBER  OR  FILE  NUMBER  WHICH  WAS  DEFINED 

C IN  SETTING  UP  THE  FILE.  ALL  OTHER  FIND  COftlENTS 

C APPLY  TO  FINDTP. 

C NOT¥ LIKE  FIND.  THE  INDEPENDENT  VARIABLE  MUST 

C BE  MONOTON ICALLY  INCREASING 

C 

C J.  ANDERSON  S-1-6S 

C EXTRAPOLATES  FOR  VALUES  OUT  OF  TABLE  RANGE 

C REMEMBERS  LAST  ARG  VALUE  AND  BEGINS  SEARCH  FROM  THAT  VALUE 

C CALLING  SEQUENCE  IS  

C 

C CALL  FINDTP  (ARG,ANS.NVAR,NARG,NANS. IDEV,NPTS. I) 

C 

C ARG  IS  THE  ARGUMENT 

C ANS  CONTAINS  RESULT  ON  EXIT 

C X IS  ONE  DIMENSIONAL  ARRAY  OF  INDEP.  VARIABLE 

C X IS  ONE  DIMENSIONAL  ARRAY  OF  DEP.  VARIABLE 

. C NPTS  IS  NUMBER  OF  TABLE  ENTRIES 

C MEM  SHOULD  BE  INITIALIZED  TO  I - AFTER  FIRST 

C . CALL  THE  SUBROUTINE  WILL  PRESERVE  VALUES  IN  IT 

C VARIABLE  WHOSE  LOCATION  IS  GIVEN  BY  NANS.  IDEV  GIVES 
C THE  DEVICE  HUMBER  OR  FILE  NUMBER  WHICH  WAS  DEFINED 

C IN  SETTING  UP  THE  FILE.  ALL  OTHER  FIND  COhf€NTB 

C APPLY  TO  FINDTP. 

C NOTE. LIKE  FIND,  THE  INDEPENDENT  VARIABLE  rUST 

C BE  MONOTON ICALLY  INCREASING 

C 

C J.  ANDERSON  6-1-65 

C EXTRAPOLATES  FOR  VALUES  OUT  OF  TABLE  RANGE 

C REMEMBERS  LAST  ARG  VALUE  AND  BEGINS  SEARCH  FROM  TRAT  VALUE 

C CALLING  SEQUENCE  ‘S  

C 

C CALL  FINDTP  (ARG,ANS,NVAR, NARG. NANS. IDEV.NPTS, 1) 

C 

C ARG  IS  THE  ARGUMENT 

C ANS  CONTAINS  RESULT  ON  EXIT 

C X IS  ONE  DIMENSIONAL  ARRAY  OF  INDEP.  VARIABLE 

C X IS  ONE  DIMENSIONAL  ARRAY  OF  PEP.  S/ARIABLE 

C NPTS  IS  NUMBER  OF  TABLE  ENTRIES 

C rCM  SHOULD  BE  INITIALIZED  TO  I - AFTER  FIRST 

C CALL  THE  SUBROUTINE  WILL  PRESERVE  VALUES  IN  IT 

C 

0002  DIMENSION  X(10).Y<10) 


ISO 


1 


0083  I READ  (IDEV'I)  (X(J).  J - UNVftR) 

0004  IF  (X(NfiRG)“flRG)  4*2,2 

0005  2 I-I-l 

0006  IF  (I-l)  3*3,1 

0007  3 I-i 
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0006  GO  TO  7 

0009  4 K • I -f  1 

0010  READ  (IDEV'K)  J - 1>NVAR) 

0OU  IF  (Y(NARG)-ARG)  5,7,7 

0012  5 I-I+l 

0013  IF  (I-NPTS)  4,6,6 

0014  6 I-NPTS- I 
tidlS  7 K • I 4-  I 

0016  READ  (IDEVM)  (X(J),  J - UNVAR) 

0017  READ  (IDEV'K)  (YCJ),  J - 1,NVAR) 

0018  ANS-X(:NANS)+(Y(NANS)-X(NANS))*(ARG-X(NARG)5/(Y(NARrO-X(NARG)) 

0019  RETURN 

0020  END 
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SUBROUTINE  FINDl  (ARG.flNS,X.Y,NPTS* I) 

C SUBROUTINE  FIND.  TABLE  LOOK-UP  UITH  LINEAR  INTERPOLATION. 

C J.  ANDERSON  6-1-65 

C EXTRAPOLATES  FOR  VALUES  OUT  OF  TABLE  RANSE 

C REMEMBERS  LAST  ARG  VALUE  AND  BEGINS  SEARCH  FROM  THAT  VALUE 

C CALL  IN  SEQUENCE  IS 

C 

C CALL  FIND  (ARG,ANS,X.Y,NPTS.rtM) 

C 

C ARG  IS  THE  ARGUMENT 

C ANS  CONTAINS  RESULT  ON  EXIT 

C X IS  ONE  DIMENSIONAL  ARRAY  OP  INDEF . VARIABLE 

C X IS  ONE  DIMENSIONAL  ARRAY  OF  DEP.  VARIABLE 

C NPTS  IS  NUMBER  OF  TABLE  ENTRIES 

C MEM  SHOULD  BE  INITIALIZED  TO  1 - AFTER  FIRST 

C CALL  THE  SUBROUTINE  DILL  PRESERVE  VALUES  IN  IT 

C 

DI^ENSION  X(10),Y(10) 
t IF  (XCn-ARG)  4.2>2 

2 I-l-l 

IF  (I-l)  3,3.1 

3 I-l 

GO  TO  7 

4 IF  (X(I+1)-ARG)  5,7,7 

5 I-I+l 

IF  <I-NPTS)  4.6,6 

6 I -NPTS- 1 

7 ANS-Y<n'KY(t+l)-Y(U)*<ARG-X{:^).’(X(«'»'l)-X(I)) 

RETURN 

END 
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SUBROUTINE  STHEQ  (A«B«D«L.h,N) 

C SUBROUTINE  SOLVING  L EQUATIONS  IN  M UNKNOUNS  WITH  N SETS  OF  RIGHT- 

C HAND  CONSTANTS.  A<L.M)  IS  THE  MATRIX  OF  COEFFICIENTS  AND  B(L,N)  IS 

C THE  MATRIX  OF  COLUMNS  OF  ANSUERS.  ON  RETURN  FROM  SUBROUTINE.  A 
C CONTAINS  QRTHOGONALIZED  COLUMNS.  B CONTAINS  THE  RESIDUALS-  AND 
C D(M.N)  CONTAINS  THE  SOLUTIONS.  FOR  MORE  EQUATIONS  THAN  UNKNOUNS 
C THE  LEAST-SQUARES  SOLUTION  IS  OBTAINED,  OTHERUISE  THE  SOLUTION  IS 
C IN  TERMS  OF  THE  FIRST  L LINEARLY  INDEPENDENT  VARIABLES.  REQUIRES.. 
C DIMENSION  ACIS.IS).  BdS.lS).  DdS.lS) 

C CALL  SIMEQ  (A.  B.  D,  L.  M.  N) 

DIMENSION  A(500,l).  B<500.1>.  C(3.3),  D(l.l) 

DO  701  I-1,M 
DO  700  J-l.M 

700  Cd.J)  - 0.0 
DO  701  K«1.N 

701  Dd.K)  • 0.0 
DO  702  J-2.M 
DO  702  I-l.L 

702  C(J.J-l)  • C(J.J-l)  + <Ad.J)i«ld.J)) 

DO  712  K-l.M 

DO  709  J-K.M 
DO  703  Il-l.L 

703  C(K.J)  - CCK.J)  + (Adl.K)iKAdl.J)) 

IF  (K-J)  706.  704.  704 

704  IF  (K-1)  709,  709,  785 

705  IF  d.E-7>*C(K.K-l)  - l.E?»C(K,K))  709.  712.  712 
786  CCK.J)  • C(K.J)/C<K.K) 

707  DO  700  12-1. L 

709  Ad2,J)  - Ad2.J)  - (Ad2.K)CCK.J)) 

709  CONTINUE 

DO  711  J2-1.N 
DO  718  13-1. L 

710  D(K.J2)  • DCK.J2)  + (A( I3,K)»0( 13. J2)/C(K.K) ) 

DO  711  14-1,.L 

711  B(I4.J2)  • fi(14.J2)  - <Ad4,K)>*iD(K.J2)) 

712  CONTINUE 

IF  (M  - 1)  715.  715,  714 

714  DO  713  1-2. M 
IT  - M+l-I 
JT  - IT+l 

DO  713  J-l.N 
DO  713  K-JT.M 

713  DCIT.J)  - DCIT.J)  - (C( IT.K)*OCK. J)) 

715  RETURN 
END 
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C CPEVES.004  22-AUG-75 
C 

0001  SUBROUTINE  EVES(N>TPRNT) 

THIS  SUBROUTINE  SOLVES  N SIMULTANEOUS  FIRST  ORDER  DIFFERENTIAL 
EQUATIONS  BY  THE  ADAMS  METHOD. 

TO  USE  ... 

- URITE  A MAIN  PROGRAM  TO  READ  THE  INPUT  DATA  AND  PERFORM  ANY 
DESIRED  AALCULATIONS  PRIOR  TO  THE  INTEGRATION.  AT  THE  POINT 
UHERE  IHE  INTEGRATION  IS  DESIRED 
CALL  EVES  (N,TPRNT) 

UHERE  N IS  THE  NUMBER  OF  DIFFERENTIAL  EQUATIONS 

TPRNT  IS  A ONE  DIMENSIONAL  ARRAY  FlLl.ED  WITH  PAIRS 
OF  VALUES  TO  CONTROL  PRINTING  (DESCRIBED  IN 
SUBROUTINE  PRINT) 

EVES  UILL  TAKE  CONTROL  AT  THIS  POINT  AND  PERFORM  THE  INTEGRATION 
UNTIL  THE  USER-PROVIDED  SUBROUTINES  INDICATE  THE  UPPER  LIMIT  HAS 
BEEN  REACHED,  (BY  SETTING  N-Q),  THEN  CONTROL  UILL  BE  RETURNED 
TO  THE  MAIN  PROGRAM 

USER  SUBROUTINES  ... 

EVES  REQUIRES  3 USER-SUBROUTINES  TO  PERFORM  THE  INTEGRATION 

-SUBROUTINE  SETUP  (T,Y.S1G,N) 

DIMENSION  T(2),Y(N).SIG(N,3) 

THIS  SUBROUTINE  SETS  THE  INITIAL  CONDITIONS  FOR  THE  INTEGRATION 
AND  DEFINES  THE  ACCURACY  REQUIRED  IN  THE  SOLUTION. 

T(l)  «■  STARTING  TIME  (ASSUMED  0.0  IF  NOT  SPECIFIED) 

T(2)  • INITIAL  TIME  INCREMENT  (ASSUMED  l.OE-5  IF  UNSPECIFIED) 
Yd)  - INITIAL  VALUES  OF  DEPENDENT  VARIABLES.  I-l.N 
(ASSUMED  0.0  IF  UNSPECIFIED) 

SIGd.l)  • REQUIRED  ACCURACY  FOR  THE  DEPENDENT  VARIABLES. 

I-l.N  UHERE  SlG-l.E-M  INDICATES  M SIGNIFICANT 
FIGURES.  (ASSUr^D  l.E-3  IF  UNSPECIFIED) 

SIG(1.2)  • MINIMUM  ABSOLUTE  ACCURACY  DESlRED-SlGd.2)^SIG(I.  I 
).  TKIS  SUSPENDS  THE  REQUIRED  NO.  OF  SIGNIFICANT 
FIGURES  IF  THE  ABSOLUTE  VALUE  OF  THE  INTEGRAL  IS 
LESS  THAN  THIS  VALUE.  ITS  USE  IS  TO  PREVENT  USE 
OF  AN  INORDINATELY  SMALL  STEP  SIZE  AS  AN 
INSETUP  (T.Y.SIG.N) 

DIMENSION  T(2).Y(N).SIG(N.3) 

THIS  SUBROUTINE  SETS  THE  INITIAL  CONDITIONS  FOR  THE  INTEGRATION 
AND  DEFINES  THE  ACCURACY  REQUIRED  IN  THE  SOLUTION. 

T(l)  - STARTING  TIME  (ASSUr^D  0.0  IF  NOT  SPECIFIED) 

T(2)  - INITIAL  TIME  UtCRECENT  (ASSUMED  l,QE-5  IF  UNSPECIFIED) 
Yd)  • INITIAL  VALUES  OF  DEPENDENT  VARIABLES.  l-l.N 
(ASSUMED  0.0  IF  UNSPECIFIED) 

SIGd.l)  - REQUIRED  ACCURACY  FOR  THE  DEPENDENT  VARIABLES. 

l-l.N  UHERE  SIG-l.E-M  INDICATES  M SIGNIFICANT 
FIGURES.  (ASSUMED  l.E-3  IF  UNSPECIFIED) 


1S5 


o n r>  o n 


SIG(I,2)  » MINIMUM  fiBSOLUTE  fiCCURftCY  DESIREO-SIG(I.2)#SIG(I.l 
).  THIS  SUSPENDS  THE  REQUIRED  NO.  OF  SIGNIFICfiNT 
FIGURES  IF  THE  ABSOLUTE  VALUE  OF  THE  INTEGRAL  IS 
LESS  THAN  THIS  VALUE.  ITS  USE  IS  TO  PREVENT  USE 
OF  AN  INORDINATELY  SMALL  STEP  SIZE  AS  AN  INTEGRAL 
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LEAVES  0.  (ASSUMED  0.  IF  UNSPECIFIED) 

SIG(I.3)  • THRESHOLD  VALUE  FOR  THE  VARIABLE  1.  EVES  UILL  HIT 
THIS  VALUE  EXACTLY  DURING  THE  INTEGRATION  AND  LET 
THE  USER  ROUTINES  KNOU  THIS  VALUE  HAS  BEEN  HIT. 
(ASSUMED  1.0E+35  IF  UNSPECIFIED) 


SUBROUTINE  DIFEQ  (T,Y.DY,N,TPR) 

DIMENSION  T(2),YCN),DY(N),TPR(2) 

THIS  SUBROUTINE  EVALUATES  THE  DERIVATIVE  VALUES  (DY(I)),AT  EACH 
PROPOSED  STEP  IN  THE  SOLUTION.  ON  EACH  CALL, TCI)  CONTAINS  THE 
CURRENT  TIME,  AND  THE  Y(I)S  ARE  THE  CURRENT  INTEGRATED  VALUES. 

THE  RESULTS  OF  A CALL  MAY  BE  REJECTED, AND  THE  STEP  SIZE  CUT  IN 
ORDER  TO  MAINTAIN  ACCURACY,  SO  NO  PERI^NENT  SUITCH  SETTING  SHOULD 
BE  MADE  IH  THIS  ROUTINE,  BASED  ON  PROPOSED  VALUES 


-SUBROUTINE  PRINT  (T, Y.DY,N,TPR) 

DIMENSION  T(2),Y(N),DY(N),TPRCN+4) 

THIS  SUBROUTINE  IS  CALLED  TO  OUTPUT  ACCEPTED  VALUES  DURING  THE 
INTEGRATION  PROCEEDURE.  THE  FREQUENCY  WITH  WHICH  IT  IS  CALLED 
IS  DEPENDENT  ON  THE  NUMBER  OF  INTEGRATION  STEPS, THE  TPRNT  ARRAY 
SUPPLIED  TO  EVES  BY  THE  MAIN  PROGRAM,  AND  THE  THRESHOLD  VALUES 
FOR  VARIABLES  SET  IN  SUBROUTINE  SETUP  AS  FOLLOWS  - 

-THE  TPRNT  ARRAY  CONSISTS  OF  PAIRS  OF  VALUES  DELTA  T -^ND  TLIM 
SUCH  THAT  DELTA  T IS  THE  PRINT  INTERVAL  UNTIL  TLIM  IS  REACHED 
, WHEREUPON  THE  NEXT  PAIR  OF  VALUES  TAKES  CONTROL. 

IF  DELTA  T .EQ.  0 EVERY  ACCEPTED  POINT  IS  PRINTED 

IF  DELTA  T .GT.  0 PRINTING  OCCURS  ONLY  AT  INTERVALS  OF 
DELTA  T 

IF  DELTA  T .LT.  8 PRINTING  OCCURS  AT  EVERY  ACCEPTED 
POINT, BUT  AMONG  THESE  POINTS  ARE 
INTERVALS  OF  ABS(DELTA  T) 

-REGARDLESS  OF  THE  TPRNT  CONTROLS, A CALL  TO  PRINT  WILL  OCCUR 
EACH  TIME  A VARIABLE  REACHES  IT'S  THRESHOLD  VALUE. 

THE  TPR  ARRAY  HAS  THE  FOLLOWING  INFORMATION  WHEN  PRINT  IS  CALLED 
TPR(l)  CURRENT  DELTA  T FROM  TPRNT  ARRAY 
TPR (2)  CURRENT  Tl  IM  FROM  TPRNT  ARRAY 
TPR(3)  0 IF  CALLED  THRU  ACCURACY  CONTROLLED  STEP  SIZE. 

+J  IF  CALLED  DUE  TO  VARIABLE  J RISING  TO  ITS 
THRESHOLD  VALUE 

-J  IF  CALLED  DUE  TO  VARIABLE  J FALLING  TO  ITS 
THRESHOLD  VALUE 

N+l  IF  CALLED  AT  SPECIFIED  PRINT  POINT 
TPR(4)  CURRENT  TIME  STEP  SIZE  FOR  INTEGRATION 

(IF  THIS  VALUE  IS  CHANGED  IN  PRINT.  THE  DIFFERENCE  TABLE 
WILL  BE  ZEROED  AND  THE  SOLUTION  RESTARTED  WITH  THE  ALTERED 
VALUE  AS  THE  STEP  SI2E.  THIS  IS  TO  GET  AROUND  EXTERNALLY 
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INDUCED  DIGCQNTINUITIES  UITHOUT  REQUIRING  THE  INTEGRATION 
TO  DISCOVER  THE  POINT  OF  OCCURRENCE  BY  HALVING  ITS  STEP.) 
TPR(5)  ROUGHEST  VARIABLE  OF  CURRENT  STEP 
TPR(6)~TPR(N-t-S)  CURRENT  THRESHOLD  SETTINGS  FOR  EACH  VARIABLE 
(ALLOUING  DYNAMIC  RESET  OF  THRESHOLDS) 
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C 

C 

0002  DIMENSION  D1  ( 12,?)  .Tl<2) , 02(12,7)  ,T2(2)  .Yl  ( 12)  ,DYl  ( 12)  >SIG(  12,3)., 

1 Y2(12).DY2(12),TPRNT(6),TPR(17) 

0003  DIMENSION  T3(2) 

0004  EQUIVALENCE  (D I ( 1)  , Y1  < 1) ) . (DU  13)  ,DYl ( 1) ) . (D2(  1) ,Y2C  I ) ) 

I , (P2(13),DY2(1)) 

0005  TQIAS"0.0 

0006  KTERR-a 

0007  I M=>N 

0008  2 DO  4 I-1,M 

0009  DO  3 J-l,7 

0010  DUI,J)-0.0 

0011  3 D2(I,J)=0.0 

0012  SIG(I.l)»1.0e-03 

0013  SIG(I,3)-1.8E+35 

0014  4 GIG(  I,2)-0.0 

0015  IPRNTM 

0016  TU1)»0.0 

0817  TU2)M.E-5 

C 

C INITIAL  SETUP  COMPLETE.  CALL  IN  THE  PROGRAl^RS  SETUP. 

C 

0910  CALL  bETUP(Tl.Yl,SlG,M) 

0019  IF  (M.LT.O)  IPRNT  - -M 

0021  .M  - N 

9022  TPR(l)-TPRNT( IPRNT) 

0023  TPR(2)»1PRNT(IPRNT+1) 

0024  •'>iEXT«T  1 ( I ) ♦ABS ( TPR  ( I ) ) 

0025  ISIJV"! 

0026  TPR (3)^8. 

0027  DO  5 I-l.« 

0028  J*I+S 

0829  3 TP»(J)rSIG(l.3) 

0830  6 CALL  DlFEQCTUYl  .DYI.M.TPR) 

0631  IF  (M)  64.  64.  7 

0032  7 CALL  PRlHT(Tl.Vl.DYl.M.TPR) 

0033  IF  (N)  64.  64.  8 

8834  8 T2(l)-Tl(l)+Tl<2) 

9035  T2(2)»TI<2) 

C 

C PREDICT  AHEAD  BY  A0AI1S  f^THOD 

C 

8836  DO  9 I-l.M 

0037  9 Y2(  I )»0. 34061  ni»D  I (1. 6)  ^O.SrS^-DUI.S) +8.41666667*0 1 (1,4)-H3.S* 

1 Dl(l.3)+Tl(?)+DY1(I)+Yl(l) 

0030  T3(1)-T2(1)+TB1AS 

0839  Ti(2)-T2(2) 

C 

C OBTAIN  TWE  DERIVATIVES  AT  THE  PROPOSED  POINT. 

C 

0040  CALL  DIFEQ(T3.Y2,DY2.M.TPR) 


1 59 


0041  701  TZSflVE— 300 

0042  IF  CM)  S4.  64,  10 

0043  10  EMAX-0. 

C 

C DIFFERENCE  ftND  COMPUTE  CRITICftL  (MAXIMUM)  ERROR  TERM. 
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0044  DO  18  I-1>M 

0045  D2(I,3)«T2(2)>k(DY2(I)-DYI(I>) 

8046  DO  U J=3.6 

0047  U D2(I,J+n-D2(I,J)-Dl(I,J) 

C 

C DELETE  DIVERGENT  DIFFERENCES. 

C 


r 

8043 

DO  15  J-3>5 

0049 

IF  (A3S(D2f.I,J+l))-flBS(D2(I,J)))  15,  12,  12 

0050 

12 

IP  <ABS(D2(I,J+2))-ttQS(D2a,J+l)))  15,  13,  13 

0051 

13 

DO  14  K"J,6 

0052 

14 

D2(LK+n“0. 

0053 

GO  TO  16 

0054 

15 

CONTINUE 

0055 

J-7 

0056 

16 

XY2  =0 . 3=«D2  ^ I , J) /(S IG  n , 1)  >i4>mXl  ( ABS  ( Y2  ( D ) , 1 . 0E-30 

0057 

IF  (flBS(W2)-EmX)  18,  17,  17 

0053 

17 

EmX=3BS(X'7Z) 

0059 

JCRflP«I 

0060 

TPRC5)=«i 

0061 

18 

CONTINUE 

L 

c 

DETERMINE  IF  ERROR  IS  UITHIN  BOUND-. 

0062 

L 

IF  (EMAX-1.0)  39,  30,  19 

L 

c 

P 

ERROR  TOO  BIG.  HALVE  THE  INTERVAL  AND  TRY  AGAIN. 

0063 

U 

p 

19 

Tl(2)“.5*Tl(2) 

L 

c 

p 

CHECK  FOR  DELTA  TIME  INSIGNIFICANT. 

0064 

L 

IF  (TKn/THZl-l.BEaG)  27,  20.  20 

3065 

p 

20 

CONTINUE 

L 

c 

•<1 

DELTA  TIME  IS  INSIGNIFICANT.  TRANSLATE  THE  ORIGIN. 

0066 

DO  21  IJA2Z»l,M 

0067 

21 

Yl(lJAZ2)-Y2(IJA2Z) 

0068 

TB I AS  “TB I AS-f-T  1 ( 1 ) ■^2 . 01KT I ( 2 ) 

0069 

T2(2)»0, 

0070 

Tl(l)-0.0 

0071 

KTERR-KTERR+l 

- 

0072 

IF  (KTERR-4)  23,  22,  22 

0073 

22 

TYPE  67,  JCRAP 

0074 

RETURN 

3075 

p 

23 

TYPE  68,  KTERR,  JCRAP 

L 

c 

p 

RESTART  THE  SOLUTION. 

0076 

L 

24 

DO  26  ;JA72-l,M 

161 


0077  DO  25  JAZ2«2,6 

0078  25  Dl(IJftZZ.Jft22)“B. 

0079  DO  26  Jfl2Z»2>7 

0080  26  D2i:iJflZZ,JftZZ)-0. 

0081  T3(l)“TBI8S+Tl(n 
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0082 

0083 

0084 

0085 

0086 
008? 
0088 

0089 

0090 


0091 

0092 

0093 

0094 

0095 

0096 

0097 

0099 

0100 
0101 
0102 

0103 

0104 

0105 

0106 

0107 

0108 
0110 

Dili 

8112 

0115 

0116 
0117 
0119 
01?.  1 
0122 
0123 
8124 
0126 
0127 

0129 

0130 

0131 


0132 

0133 


Tl(2)-l.E-5 

T3(2)-l.E-5 

CALL  DlFEQ(T3,Yl,DYl,ti.TPR) 

GO  TO  a 

27  DO  29  I-1,M 
DO  28  J-l,6 

28  D2(I,Jj^DI(I.J) 

29  D2(I,7)-0.0 
GO  TO  55 

C 

C CALCULATE  DELTA  T TO  HIT  THRESHOLD  EXACTLY 
C 

30  GO  TO  < 31,  38),ISUV 

31  IP0LD'=ABS(TPR(3J) 

TPR(3)“0. 

TRATSV-1.0 

ITERl  - 8 
DO  36  I°1,M 

IF  (I  .EQ.  IPOLD)  GO  TO  36 

IF  (SIGN(1.0, (YI(I)-TPR(I+S)))>kS1GN(1.0,(Y2(I)-TPR(I+5)))) 

1 32.-32,36 

32  TDY3=0 . •’5H<D  1 ( 1 . 6)  +0 . 33333333>t:D  1 ( 1 , 5)  +0 . 5=KD  1(  1 , 4)  +D 1 ( 1 . 3) 
TDY4-'0.,  ;66666?h<D1(1,6)+DI(I.5)+D1(I,4) 

TDY5=»1.5=KDt(I,6)+Dl(I,5) 

TDYe^’DKI.S) 

YTARG»TPR(l+5) 

!TER“0 
TROT=0.5 
XI HC -0.25 

IF  (Y2(n  -LT.  YKD)  XIKC— XINC 

763  YTEST-Y I ( I ) ■l•TRAT>i<T3C2)  >i<DY2  ( 1 1 +TDY3>vTRAT>k>k2.'2 . 0+TDY4>t‘TRAT*Y3/'6 . 0 
I -t-T DV5h<TRAT>k«4/24 . 0>«TDY6'KTR ATw5/ 1 20 . 

IF  (YTEST  -GT.  YTARG)  TRAT-TRAT-XINC 
IF  (YTEST  -LT.  YTARG)  TRAT-TPAT+XINC 
XlHC-XlNC/3. 

ITER-lTER+1 

IF  (ITER  -IT.  30)  GO  TO  700 
IF  (lTERl.EQ.0)  TRATSV  •»  VRAT 
TRATSV  “ AMINI  (TRATSV.  TROT) 

ITERl  - I 
TPR(3)«I 

IF  (Y2(n  -LT,  YKD)  TPR(3)  — I 

36  CONTINUE 

IF  (tTEPl.EQ.B)  GO  TO  39 

37  TKZ)-TSAT6V>!<Ti(2) 

151A'"2 

GO  TO  27 
C 

C TEST  FOR  ERROR  TOO  SMALL 
C 

30  T3(l)  - T2(l)  + TblAS 

IF  (EMAX  "6.0015)  39,39,48 
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c ERROR  TOO  SMALL.  ACCEPT  THIS  POINT.  BUT  DOUBLE  THE 
C INTERVAL  FOR  THE  NEXT  POINT. 

C 

0134  39  TU2)-2.0*T2(2) 
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0135 


C 

C 

C 


GO  TO  41 

ERROR  O.K.  BUY  THIS  POINT  AND  MAINTAIN  CURRENT  INTERVAL. 


0136 

0137 

0139 

0140 

0141 


40 

41 


42 


C 

C 

c 


Tli:2)«T2(2) 

IF  (TPR(3)  .NE.  0.0)  GO  TO  73 
IF  (TPR(l))  42-  44-  45 
T3<l)“T2(n+TBlAS 
T3(2)-T2(2) 


OUTPUT  THE  ACCEPTED  POINT. 


0142 

0143 

0145 

0146 

0147 

0148 


TPRC4) -T3(2) 

IF  (ABS((T3(n-TNEXT)/TNEXT) 
TZSAVE=TPR(3) 

CALL  PRINT(T3,Y2-DY2-M-TPR) 
ISW^l 

IF  (M)  64,  64,  43 


.LE.  0.00001)  TPR(3)=M+l 


C 

C 

C 


TEST  FOR  PRINTING  CONDITIONS. 


0149 

43 

IF  <T3(2)  .NE.  TPR(4))  GO  TO  200 

0151 

IF  ( ABS ( ( T3 ( 1 ) -TNEXT) /TNEXT) -0.00801) 

47.  51. 

SI 

0152 

44 

TNEXT^TPRIE) 

0153 

73 

IF  CrSd)  - TPR(2))  46.  45.  46 

0154 

45 

IF  ( ABS ( (T3 ( 1 ) -TNEXT) /TNEXT) -0 . 0000 1 ) 

99.  51, 

51 

0155 

99 

TPR(3)  = li  + 1 

0156 

TZSAVE-TPR(3) 

0157 

46 

T3(2)  - T2(2) 

0158 

T3(2)"T2(2) 

0150 

TPR(4)»r3(2) 

8160 

CALL  PRINT(T3,Y2.DY2,M.TPR) 

0161 

ISUVM 

0162 

IF  (M)  64.  64,  47 

0163 

47 

IF  (T3<2)  .NE.  TPR(4))  GO  TO  200 

0165 

IF  (T7SAVE  .EQ.  TPR(3))  TNEXT-TNEXT+A8SnPR(  1 ) ) 

0167 

IF  (TJCU-TPRIPJ-^g.  MKTPRd))  51,  48, 

48 

0168 

48 

IPRNTdPRNT+2 

0169 

TPRd)-TPRNTdPRNT) 

0170 

TPRl2)“TPRNTdPRNT+I) 

8171 

IF  dPRCD)  49.  50  . 49 

0172 

49 

TNEXT-T3d)+ABS(TPRd)) 

0173 

GO  TO  51 

0174 

50 

TNEXT-TPR(2) 

0175 

51 

IF  <Tl(2)+T3d)-TNEXT)  53.  53.  52 

0176 

52 

Tl<2)«TNEXT-T3d) 

0177 

53 

CONTINUE 

0170 

54 

T1  d)-T2d) 

0179 

55 

U"Tl(2)/T2(2) 

0180 

IF  (U-1.0)  56.  65.  56 

C 

C 


ADJUST  THE  DIFFERENCE  TABLE  FOR  INTERVAL  CHANGE. 
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0181 

0182 

0183 

0184 


56  DO  57  1-1>M 
XYZ-U-'«<U 
TMP3-D2(L3) 
ThP4-D2(I.4) 
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y, 

t 


0105  TMP5=D2(I.5) 

0186  TMP6-D2(I>6) 

0107  TMP7=D2(I,7) 

0 1 08  Tt-1F3"  (0 . 2*TriP7+0 . 25*TMP6+0 . 33333333>kTMP5+0 . 5=f<TMP4+TMP3)  >HXY2 

0189  XY2=l<Y2Hi 

0190  ThP4-  (0 , 83333333n<TMP7+0 . 9 1 6666G7>KTMP6+TMP5+TMP4)  )«XY2 

0191  XY2=XY2HJ 

0192  THP5  - ( 1 . 75=kTMP7+ 1 . 5>kTMP6+TMP5 ) HfXYZ 

0193  XYZ=LJ-'kXYZ 

0194  TMP6=*(2.0)KTMP7+TMP6)it'XYZ 

0195  T!iP7'=TMP7n<XY2*U 

0196  YKD-'i^d) 

0197  DYl(I)=nY2(I) 

0198  D 1 ( 1 , 3)  »0 . 83333333E-02:KTMP7-a . 4 1 66GG667E-0 1 )KTMP6+0 . 1 666G66G7* 
1 TfiP5-0 . 5>xTMP4+TMP3 


0199  D 1 ( I .•  4)  =0 . 58333333>KTnP6-0 . 25;t!TMP7-TMP5+0 . 5>t<TMP4+TMP3 

0200  D 1 ( I > 5)  ° 1 . 25>KTriP7- 1 . 5HiTNP6+TriP5 

0201  57  Dl(I,G)-TiiP6-2.0#TMP7 

C 

C DELETE  DIVERGENT  DIFFERENCES. 

C 

0202  58  DO  63  I-UH 

0203  DfJ  62  J-3,5 

0204  IF  (A8S(D1(I.J+1))-A0S(DI(I.J)))  G2,  59,  59 

0205  59  IF  cA8S<Dl(I.J+2))-ABS(Dl(l,JH)))  62,  60,  60 

0206  60  DO  61  K»J,5 

0207  61  Dl(I.K+l)-0.0 

0208  GO  TO  63 

0209  62  CONTINUE 

0210  63  CONTINUE 

0211  GO  TO  8 

0212  64  RETURN 

0213  S5  Tl(l)-T2(n 

0214  Tl(2)-T2(2) 

0215  DO  66  I-l,M 

0216  DO  66  J-1.6 

0217  66  Dl< I,J)-D2(I,J) 

0218  GO  TO  8 

0219  67  FORMAT  ( lH1.30X.50Htt)RE  THAN  THREE  ORIGIN  TRANSLATION  DUE  TO  VAHIA 

IDLE. 13./5X. 12HRUH  ABORTED.) 

0220  C8  FORMAT  ( IH0/'1M0.30X.22HORIGIN  TRANSLATION  NO. , I3.26X, 15HCUE  TO  VAR 

IIABLE.I3./IH8) 

C 

C CLEAR  DIFFERENCE  TABLE  AT  CUSTOMERS  REQUEST 

C 

0221  200  DO  201  IJA22-I.M 

0222  DO  202  JA22»2.6 

0223  202  Dt( IJA22,JA23)-0.0 

0224  DO  201  JAZZ«2,? 

0225  201  D3(!JA22.JA22)-0.a 

0226  Tl{i)-T2(l) 

022?  DO  204  I-l.M 
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0:|28  204  Yl(I)-Y2(I) 

0229  Tl(2)-TPR(4) 

0230  CALL  DIFEQ  (T3,Yl,DYl>ti,TPR) 

0231  GO  TO  0 

0232  END  . 


! 
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ATTN : Herman  Krier 
Transporation  Bldg,  Rm  105 
Urbana,  IL  61801 

1 University  of  Massachusetts 
Dept  of  Mechanical  Engineering 
ATTN:  Karl  Jakus 
Amherst.  MA  01002 

1 University  of  Utah 

Dept  of  Chemical  Engineering 
ATTN;  Alva  D.  Baer 
Park  Eldg.,  Room  307 
Salt  Lake  City,  UT  84112 

1 University  of  Waterloo 

Dept  of  Mechanical  Engineering 
ATTN:  Clarke  E.  Herraance 
Waterloo,  Ontario,  Canada 


Aberdeen  Proving  Ground 

Marine  Corps  Ln  Ofc 
Dir,  USAMSAA 
Dir,  USAMTD 

AITN:  Mr.  H.  Anderson 

Mr.  W.  Rieden 

Mr.  S.  Walton 

Cdr , USATF.COM 

ATTN;  DRSTE-FA,  Mr.  J.  Byrne 


1 University  of  Delaware 
A'n’N:  T.  C.  Brill 
Dep.artment  of  Chemistry 
Newark,  DE  19711 
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